!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This file contains all the algorithms that are used
! to solve for the stationary general equilibrium of the
! economy, to simulate a panel of firms, to construct
! the model counterparts of targeted moments and other
! quantities of interest (e.g., the levels of output,
! consumption, etc.), and to conduct quantitative 
! decompositions.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !sa call the subroutine fcn
    !fcn is the program which generate loss function once given parameter values.

    subroutine fcn(nprmt,params,func,poe)

    use parameters
    use globals
    use tools
    !use svd
    use datatype




    INTEGER, INTENT(IN)  ::poe,nprmt
    real*8,intent(inout)::params(nprmt)
    real*8,allocatable:: datamom(:,:),simmom(:,:),tsimmom(:,:),weight(:,:),modelmom2(:,:),modelmom2_w(:,:),&
        ww(:,:),datamom_w(:,:),simmom_w(:,:),tsimmom_w(:,:),weight_w(:,:),aggmom(:,:),aggmom_w(:,:),pickout(:,:)
    real*8,intent(out):: func
    integer::errorflag


    integer ::    nfirm, ngv, cluster, samplesize
    real(8) :: ppp(nprmt),  gwg(nprmt,nprmt), wsw(nmom,nmom), gwswg(nprmt,nprmt), igwg(nprmt,nprmt), vc(nprmt,nprmt), standarderror(nprmt), gigwgg(nmom,nmom), &
        gw(nmom,nprmt), eye(nmom,nmom), eyegg(nmom,nmom), vpe(nmom,nmom), ivpe(nmom,nmom), nsim, testgv, jtest, outtest
    character (len = 50) :: gvfil
    character (len = 50) :: momdfil
    character (len = 50) :: resfil
    real(8), allocatable :: w(:,:), momdiff(:), allmoms(:), simmoms(:), datamoms(:), swi(:,:), sw(:,:), gv(:), gvi(:),tstat(:), &
        eyew(:,:), gvall(:), gr(:,:), outg(:,:), mominfl(:,:), covmx(:,:)

    character (len=28) :: covfil

    allocate (w(nmom,nmom))
    allocate (swi(nmom,nmom))
    allocate (sw(nmom,nmom))
    allocate (eyew(nmom,nmom))
    allocate (simmoms(nmom))
    allocate (gr(nmom,nprmt))
    allocate (momdiff(nmom))
    allocate (datamoms(nmom))
    allocate (tstat(nmom))

    allocate(datamom(nmom,1))
    allocate(modelmom2(nmom,1))
    allocate(tsimmom(nmom,1))
    allocate(simmom(nmom,1))
    allocate(weight(nmom,nmom))
    allocate(ww(nmom,nmom))
    allocate(pickout(nmom,1))

    w = 0d0
    sw = 0d0
    swi = 0d0
    datamoms = 0d0
    simmoms = 0d0
    momdiff = 0d0





    if (statistics==0) then
        call policy(params) ! given parameters, solve the policy rule
        call makemoments(modelmom2)    !calculate model moments
    end if


    !    !============= Read in the data moments ===================
    !      allmoms = 0.0
    !      momdfil = "Datfils/mom.asc"
    !      open (unit=1253,file=trim(momdfil),status="old")
    !      do jj=1,allnmom
    !       read(1253,"(f26.16)") allmoms(jj)
    !      enddo
    !      close (unit=1253)
    !
    !      datamoms = allmoms(pickmom)



    if (setting==1) then
        datamom(1,1)=.4717088 !mean equity ratio
        datamom(2,1)=0.0243245000 !beta1(equity,innov_stock)
        datamom(3,1)=0.0060232000 !beta1(merger,innov_stock)
        datamom(4,1)=-0.0006818000 !beta2(merger,innov_stock)
        datamom(5,1)=0.0191241000 !avg_merger_prob
        datamom(6,1)=-0.0295337000 !beta1(completion,innov_stock)
        datamom(7,1)=0.0316473000 !avg_realized_gain
        datamom(8,1)=0.0108856000 !avg_value_loss
        datamom(9,1)=0.337  !avg relative size
        datamom(10,1)=0.0849727000 !avg r&d intensity
        datamom(11,1)=0.9443351 !rho(log(mve))
        datamom(12,1)=0.36413854!cv (log_mve )
        datamom(13,1)=0.0520995 !firm growth rate
        datamom(14,1)=0.0457087 !firm entry rate


    end if


    pickout(:,1)=1d0


    simmom=(1d0-modelmom2/datamom)*pickout
    tsimmom=simmom



    if (complicated == 1) then


        cluster = 0
        call readweightmatrix(cluster,w) !read in weighting matrix

        cluster = 1  !read in the covariance matrix omega
        call readweightmatrix(cluster,sw)


        !    if (cluster == 1) then
        !        covfil = "Datfils/omega.txt"
        !        open (unit=62,file=trim(covfil),status="old")
        !    else
        !       ! covfil = "Datfils/covmx.txt"
        !        covfil = "Datfils/Wcovmx.txt"
        !        open (unit=62,file=trim(covfil),status="old")
        !    endif



    elseif (complicated == 0) then


        cluster = 0
        call readweightmatrix(cluster,w)
        call findinv(w, sw, nmom, errcode)


        !            cluster = 1
        !            call readweightmatrix(cluster,sw)
        !            call findinv(sw, w, nmom, errcode)

    endif


    
    weight=0d0
    
    do i=1,nmom
        weight(i,i)=1d0
    end do

    simmom=(1d0-modelmom2/datamom)*pickout

    tsimmom=simmom

    simmom=matmul(weight,simmom)

    !given data moments and model moments, generate loss function as output

    func=sum(tsimmom*simmom)


    if (diffv>tolerv.OR.diff_dist > toler_dist.or. diff_pol>0.1) then
        func = 10000000000.0
    endif


    if (report_moments==1) then
        print *, 'value of this iteration is', func

        if (estimate==0) then
            OPEN (UNIT=1,FILE='datamom.txt',STATUS='replace')
            DO i=1,nmom
                WRITE(1,FMT='(f20.15)') datamom(i,1)
            END DO
            CLOSE (UNIT=1)

            print 110, 'Moments to match','        datamom','       modelmom','      pickout'
            print 111, 'mean_equity_ratio    ',datamom(1,1),modelmom2(1,1),pickout(1,1)
            print 111, 'beta1(equity,zt)     ',datamom(2,1),modelmom2(2,1),pickout(2,1)
            print 111, 'beta1(merger,zt)     ',datamom(3,1),modelmom2(3,1),pickout(3,1)
            print 111, 'beta2(merger,zt)     ',datamom(4,1),modelmom2(4,1),pickout(4,1)
            print 111, 'avg_merger_prob      ',datamom(5,1),modelmom2(5,1),pickout(5,1)
            print 111, 'beta1(completion,zt) ',datamom(6,1),modelmom2(6,1),pickout(6,1)
            print 111, 'avg_realized_gain    ',datamom(7,1),modelmom2(7,1),pickout(7,1)
            print 111, 'avg_value_loss       ',datamom(8,1),modelmom2(8,1),pickout(8,1)
            print 111, 'avg_rel_size         ',datamom(9,1),modelmom2(9,1),pickout(9,1)
            print 111, 'avg_rnd_intensity    ',datamom(10,1),modelmom2(10,1),pickout(10,1)
            print 111, 'rho(log(mve))        ',datamom(11,1),modelmom2(11,1),pickout(11,1)
            print 111, 'cv(log(mve))         ',datamom(12,1),modelmom2(12,1),pickout(12,1)
            print 111, 'firm growth rate     ',datamom(13,1),modelmom2(13,1),pickout(13,1)
            print 111, 'firm entry rate      ',datamom(14,1),modelmom2(14,1),pickout(14,1)



110         format (4A)
111         format (A,f10.4,f10.4,f10.1)
112         format (A,f10.4,f10.4)



            OPEN (UNIT=1000,FILE='writeout_xls.txt',STATUS='replace')



            do i=1,nmom
                write(1000,*) modelmom2(i,1)
            end do
            write(1000,*), output
            write(1000,*), bbeta8(2,1)
            write(1000,*), avg_innovation
            write(1000,*), avg_i_policy
            write(1000,*), avg_growth_rate
            write(1000,*)
            WRITE(1000,*)
            write(1000,*), bbeta6(2,1)
            write(1000,*), avg_value_loss
            write(1000,*), avg_realized_gain_dollar
            write(1000,*), avg_realized_gain_dollar_unconditional


            write(1000,*)
            do i=1,nparam
                WRITE(1000,*) params(i)
            end do
            write(1000,*), consumption


            close(1000)

        else


            print *, 'datamom'
            print *, datamom
            print *, 'modelmom'
            print *, modelmom2
            print *, 'parameters are'
            print *, params
            !100 format (f10.4)
        end if
    end if



    modelmom(1:nmom,1)=modelmom2(:,1)
    modelmom_nvary(:)=modelmom2(:,1)




    deallocate (w)
    deallocate (swi)
    deallocate (sw)
    deallocate (eyew)
    deallocate (simmoms)
    deallocate (gr)
    deallocate (momdiff)
    deallocate (datamoms)
    deallocate (tstat)

    deallocate(datamom)
    deallocate(modelmom2)
    deallocate(tsimmom)
    deallocate(simmom)
    deallocate(weight)
    deallocate(ww)
    deallocate(pickout)



    end subroutine fcn












    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!POLICY!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


    subroutine policy(params)

    USE parameters
    USE globals
    use tools
    use omp_lib

    IMPLICIT NONE


    !include "mkl_blas.fi"
    !include "mkl_lapack.fi"

    integer*8:: int1, int2
    real*8,intent(in):: params(nparam)
    real*8          ::  cmid,fmin,fmax

    real*8:: cv,outside_value
    integer::s_idx,i_idx,ii_idx,j_idx
    real*8::g1s0,g2s0
    integer ( kind = 4 ) statusvalue
    integer ( kind = 4 ) whichvalue
    real*8::Gvalue,  pvalue, qvalue,xvalue,boundvalue,pvalue2, qvalue2,xvalue2
    real*8:: meanvalue, sdvalue
    integer::errorflag



    alpha=params(1)
    tau=1d0-params(1)
    gamma1=params(2)
    info_asym=params(3)
    LAMBDA1=params(4)
    chi=params(5)
    stdea2=params(6)
    rhoa2=params(7)
    mu_ub=params(8)
    psi=params(9)





    xlow=0d0
    eta=1d0
    stdea= 0.55228
    rhoa=0.84163
    phi_c=2d0

    if (case_type==1) then
        omega=0d0
        omega2=0d0
    else
        omega=0d0
        omega2=0d0
    end if

    if (case_type==2) then
        outside_option=0.01
    elseif (case_type==3) then
        outside_option=0.75
    else
        outside_option=1d0
    end if

    gammae2=1d0


 
    F_exog = 0d0
    F_exog(1,1) = 1d0


    pa=0.5
    pt=0.5
    c=0.001 
    s=0d0 
    cmin=0.00001

    if (illustration==0) then
        cmax=10000d0 
    else
        cmax=100d0 
    end if

    eps=0.2
    iter_dist = 0
    diffmu=1d0
    diff_dist = 1000d0





    if (even==1) then
        do i=1,numc
            c(i)=cmin+(cmax-cmin)/(real(numc,8)-1d0)*(real(i,8)-1d0)
        end do
        do i=1,nums
            s(i)=smin+(smax-smin)/(real(nums,8)-1d0)*(real(i,8)-1d0)
        end do
    elseif (even==2) then

        do i=1,numc
            c(i)=exp((log(cmax)-log(cmin))/(real(numc,8)-1d0)*(real(i,8)-1d0)+log(cmin))
        end do

   
        do i=1,nums
            if (i<nums-scut) then
                s(i)=exp((log(smax)-log(smin))/(real(nums,8)-1d0)*(real(i,8)-1d0)+log(smin))
            else
                s(i)=s(nums-scut-1)+(smax-s(nums-scut-1))/(scut+2-1d0)*(real(i-nums+scut+2,8)-1d0)

            end if

        

        end do



    end if


    s1(1,:)=s(:)
    c1(1,:)=c(:)



    x(1)=xlow



    muza=0d0
    call markov(rhoa,stdea**2d0,numz,str,muza,muza,transz,za)



    do i=1,numi
        mu_value(i)=(i-1)*mu_ub*2d0/(numi-1)   -mu_ub/2d0
    end do

    i_cost=chi*(mu_value-mu_value(1))**phi_c


    do i=1,numi
        muza=mu_value(i)
        call markov_feed(rhoa2,stdea2**2d0,numz,za,muza,muza,transz)


        p2(3,1) = 0d0
        p2(2,1) = transz(1,1)
        p2(1,1) = transz(1,2)
        p2(:,1) = p2(:,1)/sum(p2(:,1))

        do ii=2,numz-1
            p2(3,ii) = transz(ii,ii-1)
            p2(2,ii) = transz(ii,ii)
            p2(1,ii) = transz(ii,ii+1)
            p2(:,ii) = p2(:,ii)/sum(p2(:,ii) )
        enddo

        p2(3,numz) = transz(numz,numz-1)
        p2(2,numz) = transz(numz,numz)
        p2(1,numz) = 0d0
        p2(:,numz) = p2(:,numz)/sum(p2(:,numz))

        p3(:,:,i)=p2(:,:)

    end do



    za=exp(za)



    c_threshold=0d0
    error_c=0d0
    c_threshold2=0d0
    error_c2=0d0
    gc=0d0
    gs=0d0
    gc_old=0d0
    gs_old=0d0



    V=1d0 
    Vnew=1d0 
    Va=0d0 
    Vanew=0d0 
    Vt=0d0 
    Vtnew=0d0 
    Vtnew_h=0d0 
    Vtnew_m=0d0 
    Vtnew_l=0d0 
    V_interp=1d0 
    be_acq=0d0 
    V=1d0 
    Vnew=1d0 
    Va=0d0 
    Vanew=0d0 
    Vt=0d0 
    Vtnew=0d0 
    Vtnew_h=0d0 
    Vtnew_m=0d0 
    Vtnew_l=0d0 
    V_interp=1d0 
    be_acq=0d0 
    Fs =1d0/(numz*numx) 

    v_prod = 0d0
    v_prod_new = 0d0
    i_policy=1
    i_policy_new=1
    Vanew_i=0d0
    Vtnew_h_i=0d0 
    Vtnew_m_i=0d0 
    Vtnew_l_i=0d0 
    V_meet_i=0d0
    Vnew_i=0d0
    indexc_i=1
    indexs_i=1
    offermade_i = 1
    gc_i = 0d0
    gs_i = 0d0




    Fa = 0d0 
    Ft = 0d0 
    Ftp = 0d0

    f=0d0 

    Va_2=0d0 
    indexs=0d0 
    indexc=0d0 
    indexs_old=0d0 
    indexc_old=0d0 
    gacp=0d0 
    gacp2=0d0 
    vv_record2=0d0 
    gacp_record=0d0 
    temp1new_record=0d0 
    offermade=0d0 
    offermade_old=0d0 
    gacp_h=0d0 
    gacp_m=0d0 
    gacp_l=0d0 

    theta=0.5

    ! !$OMP PARALLEL PRIVATE(i,j)
    ! !$OMP DO SCHEDULE(dynamic)
    do i = 1,numz
        do j = 1,numx
            if (i <= (numz/2))then
                be_acq_guess(i,j) = 0
            else
                be_acq_guess(i,j) = 1
            end if
        end do
    end do
    !!$omp end  do
    !!$omp end parallel

    !   !$OMP PARALLEL PRIVATE(i,j,k,ip,iip)
    !  !$OMP DO SCHEDULE(dynamic)
    do i = 1,numz
        do j = 1,numx
            if (be_acq_guess(i,j) == 1 ) then! Acquirer
                do k = -1,1 ! acquirer innovation realization
                    ip = max(min(i+k,numz),1)
                    Fa(ip,j) = Fa(ip,j) + p2(2-k,i)*(Fs(i,j)/sum(Fs))
                end do
            else ! Target
                Ft(i,j) = Ft(i,j) + (Fs(i,j)/sum(Fs))
                do k = -1,1 ! target innovation realization
                    iip = max(min(i+k,numz),1)
                    Ftp(iip,j) = Ftp(iip,j) + p2(2-k,i)*(Fs(i,j)/sum(Fs))
                end do
            end if
        end do
    end do

    !  !$omp end  do
    !  !$omp end parallel




    theta=sum(Fa)/sum(Ft)

    pa=min(min(eta*theta**(mu-1d0),1d0),1d0/theta)
    pt=min(min(eta*theta**mu,1d0),theta)




    zt=za

    za1(:,1)=za



    do i=1,numz
        do j=1,numz
            f(i,j)=gamma1*za(i)**alpha*zt(j)**tau
        end do
    end do


    do i=1,numz
        do j=1,numz
            gap =minval(abs(f(i,j)-za1(:,1)),1)
            ttt =minloc(abs(f(i,j)-za1(:,1)),1)
            if (f(i,j)>za(ttt)) then
                i_lb=ttt
                i_ub=min(ttt+1,numz)
                if (i_ub>i_lb) then
                    p_lb=1d0-gap/(za(i_ub)-za(i_lb))
                    p_ub=1d0-p_lb
                else
                    p_lb=1d0
                    p_ub=0d0
                end if
            else
                i_lb=max(ttt-1,1)
                i_ub=ttt
                if (i_ub>i_lb) then
                    p_lb=gap/(za(i_ub)-za(i_lb))
                    p_ub=1d0-p_lb
                else
                    p_lb=0d0
                    p_ub=1d0
                end if
            end if
            f2(i,j)=p_lb*za(i_lb)+p_ub*za(i_ub)

        end do
    end do




    f=f2


    
    maxiter2=500

    allocate(V_record(numz,numx,maxiter2))

    V_record=0d0
    ! Fs_record_all=0d0


    zbar = 0d0
    do i=1,numz
        do j=1,numx
            zbar = zbar + zt(i)*Fs(i,j)
        enddo
    enddo



    diffv=1d0


    !covfil = "i_policy_fixed.txt"
    !
    !    do ii = 1, numi
    !        read(62,*) aaa(:, 1)
    !    enddo

    if (GE_fixed_Fs_innov==1) then

        covfil2 = "Datfils/i_policy_fixed.txt"
        open (unit=62,file=trim(covfil2),status="old")
        do ii = 1, numz
            read(62,*) ipol_load(ii, 1)
        enddo


        i_policy=ipol_load
        i_policy_new=ipol_load

        covfil2 = "Datfils/Fs_fixed.txt"
        open (unit=62,file=trim(covfil2),status="old")
        do ii = 1, numz
            read(62,*) Fs_load(ii, 1)
        enddo



        Fs=Fs_load
        Fs_old=Fs_load
    end if




    do while (iter_dist <= max_iter_dist .and. (diff_dist > toler_dist .or. diffv>tolerv ))




        if (diffv>1e-4) then
            maxiter2=500 
        else
            maxiter2=500 
        end if
        diffv=1d0
        counter=1

        pi=0.66490574  ! A*zeta*(kappa/(1d0/beta-1d0+delta))**(kappa/(1d0-kappa))
        k_helper=(pi/zeta)**(1d0/kappa)

        do while (diffv>tolerv.and.counter<maxiter2)
            !CALL CPU_TIME(t1a)
            !call system_clock(int1)

            !$OMP PARALLEL PRIVATE(i,j,jj,gap,ttt,i_lb,i_ub,p_lb,p_ub)
            !$OMP DO collapse(3) SCHEDULE(dynamic)
            do i=1,numz
                do j=1,numz
                    do jj=1,numx


                        gap =minval(abs(f(i,j)-za1(:,1)),1)
                        ttt =minloc(abs(f(i,j)-za1(:,1)),1)
                        if (f(i,j)>za(ttt)) then
                            i_lb=ttt
                            i_ub=min(ttt+1,numz)
                            if (i_ub>i_lb) then
                                p_lb=1d0-gap/(za(i_ub)-za(i_lb))
                                p_ub=1d0-p_lb
                            else
                                p_lb=1d0
                                p_ub=0d0
                            end if
                        else
                            i_lb=max(ttt-1,1)
                            i_ub=ttt
                            if (i_ub>i_lb) then
                                p_lb=gap/(za(i_ub)-za(i_lb))
                                p_ub=1d0-p_lb
                            else
                                p_lb=0d0
                                p_ub=1d0
                            end if
                        end if
                        V_interp(i,j,jj)=V(i_lb,jj)*p_lb+V(i_ub,jj)*p_ub

                        !                        if (f(i,j)>za(numz)) then
                        !                            V_interp(i,j,jj)=V(numz-1,jj)+(V(numz,jj)-V(numz-1,jj))*(f(i,j)-za(numz-1))/(za(numz)-za(numz-1))
                        !                        end if

                    end do
                end do
            end do
            !$omp end  do
            !$omp end parallel
            !CALL CPU_TIME(t1b)

            !print *, 't1b-t1a',t1b-t1a

            !call system_clock(int2)

            !print *, 'First Part',(int2-int1)

            !call system_clock(int1)
			
			! Raw cash values in the code correspond to c*(1-s) in the draft due to legacy model specification.

            !$OMP PARALLEL PRIVATE(i,j,ii,jj,is,temp2s,temp2new,diff,iter1,cmid,cmax, cmin,fmin,fmax,cv,outside_value,s_idx,i_idx,ii_idx,j_idx,g1s0,g2s0,temp1h,temp1l,temp1_record,x_vp)
            !$OMP DO collapse(5) SCHEDULE(dynamic)
            do i=1,numz  !acquirer
                do j=1,numx
                    do ii=1,numz !target
                        do jj=1,numx

                            do is=1,nums

                                temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                temp2s=max(temp2s,1e-10)
                                temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                cmin=0d0 !-0.1
                                cmax=c(numc)+0.1
                                diff=1d0
                                iter1=1


                                do while (abs(diff)>toler.and.iter1<maxiter)
                                    cmid=(cmin+cmax)/2d0


                                    cv = cmid
                                    outside_value = temp2new
                                    s_idx = is
                                    i_idx = i
                                    ii_idx = ii
                                    j_idx = j

                                    x_vp = ((zt(ii_idx)/zt(i_idx))**phi * (zt(ii_idx)+zt(i_idx))**phi2)
                                    g1s0=LAMBDA1*( 1d0-cv/(cv+s(s_idx)*(  max(0d0,pi*f(i_idx,ii_idx)-x(j_idx)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i_idx,ii_idx,j_idx)*fv)      )) )**gammae2*f(i_idx,ii_idx)
                                    g1s0=max(g1s0,0d0)
                                    g2s0=LAMBDA1*( 1d0-cv/(cv+s(s_idx)*(  max(0d0,pi*f(max(i_idx-1,1),ii_idx)-x(j_idx)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i_idx-1,1),ii_idx,j_idx)*fv)     )) )**gammae2*f(max(i_idx-1,1),ii_idx)
                                    g2s0=max(g2s0,0d0)
                                    temp1h = cv+s(s_idx)*(pi*f(i_idx,ii_idx)-x(j_idx)*x_vp-g1s0+(1d0-psi)/(1d0+r)*V_interp(i_idx,ii_idx,j_idx))
                                    temp1l = cv+s(s_idx)*(pi*f(max(i_idx-1,1),ii_idx)-x(j_idx)*x_vp-g2s0+(1d0-psi)/(1d0+r)*V_interp(max(i_idx-1,1),ii_idx,j_idx))
                                    temp1h = max(temp1h,1e-10)
                                    temp1l = max(temp1l,1e-10)
                                    temp1_record = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)
                                    diff = temp1_record - outside_value



                                    cv = cmin
                                    outside_value = temp2new
                                    s_idx = is
                                    i_idx = i
                                    ii_idx = ii
                                    j_idx = j

                                    x_vp = ((zt(ii_idx)/zt(i_idx))**phi * (zt(ii_idx)+zt(i_idx))**phi2)
                                    g1s0=LAMBDA1*( 1d0-cv/(cv+s(s_idx)*(  max(0d0,pi*f(i_idx,ii_idx)-x(j_idx)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i_idx,ii_idx,j_idx)*fv)      )) )**gammae2*f(i_idx,ii_idx)
                                    g1s0=max(g1s0,0d0)
                                    g2s0=LAMBDA1*( 1d0-cv/(cv+s(s_idx)*(  max(0d0,pi*f(max(i_idx-1,1),ii_idx)-x(j_idx)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i_idx-1,1),ii_idx,j_idx)*fv)     )) )**gammae2*f(max(i_idx-1,1),ii_idx)
                                    g2s0=max(g2s0,0d0)
                                    temp1h = cv+s(s_idx)*(pi*f(i_idx,ii_idx)-x(j_idx)*x_vp-g1s0+(1d0-psi)/(1d0+r)*V_interp(i_idx,ii_idx,j_idx))
                                    temp1l = cv+s(s_idx)*(pi*f(max(i_idx-1,1),ii_idx)-x(j_idx)*x_vp-g2s0+(1d0-psi)/(1d0+r)*V_interp(max(i_idx-1,1),ii_idx,j_idx))
                                    temp1h = max(temp1h,1e-10)
                                    temp1l = max(temp1l,1e-10)
                                    temp1_record = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)
                                    fmin = temp1_record - outside_value




                                    cv = cmax
                                    outside_value = temp2new
                                    s_idx = is
                                    i_idx = i
                                    ii_idx = ii
                                    j_idx = j

                                    x_vp = ((zt(ii_idx)/zt(i_idx))**phi * (zt(ii_idx)+zt(i_idx))**phi2)
                                    g1s0=LAMBDA1*( 1d0-cv/(cv+s(s_idx)*(  max(0d0,pi*f(i_idx,ii_idx)-x(j_idx)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i_idx,ii_idx,j_idx)*fv)      )) )**gammae2*f(i_idx,ii_idx)
                                    g1s0=max(g1s0,0d0)
                                    g2s0=LAMBDA1*( 1d0-cv/(cv+s(s_idx)*(  max(0d0,pi*f(max(i_idx-1,1),ii_idx)-x(j_idx)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i_idx-1,1),ii_idx,j_idx)*fv)     )) )**gammae2*f(max(i_idx-1,1),ii_idx)
                                    g2s0=max(g2s0,0d0)
                                    temp1h = cv+s(s_idx)*(pi*f(i_idx,ii_idx)-x(j_idx)*x_vp-g1s0+(1d0-psi)/(1d0+r)*V_interp(i_idx,ii_idx,j_idx))
                                    temp1l = cv+s(s_idx)*(pi*f(max(i_idx-1,1),ii_idx)-x(j_idx)*x_vp-g2s0+(1d0-psi)/(1d0+r)*V_interp(max(i_idx-1,1),ii_idx,j_idx))
                                    temp1h = max(temp1h,1e-10)
                                    temp1l = max(temp1l,1e-10)
                                    temp1_record = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)
                                    fmax = temp1_record - outside_value


                                    if (fmax>=fmin) then
                                        if (diff>=0d0) then
                                            cmax=cmid
                                        else
                                            cmin=cmid
                                        end if
                                    else


                                        print *, 'XXXXXXX'
                                        print *, cmin,cmid,cmax,s(is)
                                        print *, fmin,diff,fmax
                                        if (diff>0d0) then
                                            cmin=cmid
                                        else
                                            cmax=cmid
                                        end if
                                    end if
                                    iter1=iter1+1

                                end do

                                if (diff > -1e-2) then
                                    c_threshold2(i,j,ii,jj,is)=cmid
                                    error_c2(i,j,ii,jj,is)=diff
                                else

                                    c_threshold2(i,j,ii,jj,is)=1e8
                                    error_c2(i,j,ii,jj,is)=diff
                                end if


                            end do
                            !   !$omp end  do
                            !  !$omp end parallel
                        end do
                    end do
                end do
            end do

            !$omp end  do
            !$omp end parallel
            !call system_clock(int2)

            !print *, 'Second Part',(int2-int1)

            !call system_clock(int1)

            !acquirer problem

            !$OMP PARALLEL PRIVATE(i,j,ii,jj,inn,iii,ss,cc,cc1,cc2,cc3,g1,g2,temp1,temp2,temp1b,temp2b,vv_h,vv_m,vv_l,vv,vvc,vvs,idxc,idxs,p2,x_vp,offermade,Vanew, gc,gs,indexc,indexs)
            !$OMP DO collapse(4) SCHEDULE(dynamic)
            do i=1,numz  !acquirer
                do j=1,numx !acquirer
                    do ii=1,numz !observed zt_hat
                        do jj=1,numx

                            ss=reshape(s1,(/3,nums/),s1,(/2,1/))


                            cc=0d0

                            cc(1,:) = c_threshold2(i,j,min(ii+1,numz),jj,:)
                            cc(2,:) = c_threshold2(i,j,ii,jj,:)
                            cc(3,:) = c_threshold2(i,j,max(ii-1,1),jj,:)
                            where (cc<0d0)
                                cc=0d0
                            end where


                            cc1(1,:) = cc(1,:)
                            cc1(2,:) = cc(1,:)
                            cc1(3,:) = cc(1,:)

                            cc2(1,:) = cc(2,:)
                            cc2(2,:) = cc(2,:)
                            cc2(3,:) = cc(2,:)

                            cc3(1,:) = cc(3,:)
                            cc3(2,:) = cc(3,:)
                            cc3(3,:) = cc(3,:)

                            where (cc>=cc1)
                                gacp_h(i,j,ii,jj,:,:)=1d0  !check this
                            elsewhere
                                gacp_h(i,j,ii,jj,:,:)=0d0  !check this
                            end where

                            where (cc>=cc2)
                                gacp_m(i,j,ii,jj,:,:)=1d0
                            elsewhere
                                gacp_m(i,j,ii,jj,:,:)=0d0  !check this

                            end where

                            where (cc>=cc3)
                                gacp_l(i,j,ii,jj,:,:)=1d0
                            elsewhere
                                gacp_l(i,j,ii,jj,:,:)=0d0  !check this

                            end where

                            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                            ! Imperfect Information Case
                            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                            iii = min(ii+1,numz)
                            x_vp = ((zt(iii)/zt(i))**phi * (zt(iii)+zt(i))**phi2)
                            g1=LAMBDA1*( 1d0-cc/(cc+ss*(  max(0d0,pi*f(i,iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)*fv)  )     ) )**gammae2*f(i,iii) 

                            where (g1<0d0)
                                g1=0d0
                            end where



                            g2=LAMBDA1*( 1d0-cc/(cc+ss*(   max(0d0,pi*f(max(i-1,1),iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)*fv)  )       ) )**gammae2*f(max(i-1,1),iii) 
                            where (g2<0d0)
                                g2=0d0
                            end where

                            temp1=pi*f(i,iii)-x(j)*x_vp-g1
                            temp2=pi*f(max(i-1,1),iii)-x(j)*x_vp-g2
                            temp1b=(-cc+(1d0-ss)*(temp1+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)))
                            temp2b=(-cc+(1d0-ss)*(temp2+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)))


                            where (temp1b<1e-10)
                                temp1b=1e-10
                            end where
                            where (temp2b<1e-10)
                                temp2b=1e-10
                            end where


                            vv_h=gacp_h(i,j,ii,jj,:,:)*(kesi*temp1b(:,:)**(1d0-omega)/(1d0-omega)+(1d0-kesi)*temp2b(:,:)**(1d0-omega)/(1d0-omega))&
                                +(1d0- gacp_h(i,j,ii,jj,:,:))*((pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))




                            iii=ii
                            x_vp = ((zt(iii)/zt(i))**phi * (zt(iii)+zt(i))**phi2)
                            g1=LAMBDA1*( 1d0-cc/(cc+ss*(  max(0d0,pi*f(i,iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)*fv)  )     ) )**gammae2*f(i,iii) 
                            where (g1<0d0)
                                g1=0d0
                            end where
                            g2=LAMBDA1*( 1d0-cc/(cc+ss*(   max(0d0,pi*f(max(i-1,1),iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)*fv)  )       ) )**gammae2*f(max(i-1,1),iii) 
                            where (g2<0d0)
                                g2=0d0
                            end where
                            temp1=pi*f(i,iii)-x(j)*x_vp-g1
                            temp2=pi*f(max(i-1,1),iii)-x(j)*x_vp-g2
                            temp1b=(-cc+(1d0-ss)*(temp1+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)))
                            temp2b=(-cc+(1d0-ss)*(temp2+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)))

                            where (temp1b<1e-10)
                                temp1b=1e-10
                            end where
                            where (temp2b<1e-10)
                                temp2b=1e-10
                            end where

                            vv_m=gacp_m(i,j,ii,jj,:,:)*(kesi*temp1b(:,:)**(1d0-omega)/(1d0-omega)+(1d0-kesi)*temp2b(:,:)**(1d0-omega)/(1d0-omega))&
                                +(1d0- gacp_m(i,j,ii,jj,:,:))*((pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))



                            iii = max(ii-1,1)
                            x_vp = ((zt(iii)/zt(i))**phi * (zt(iii)+zt(i))**phi2)
                            g1=LAMBDA1*( 1d0-cc/(cc+ss*(  max(0d0,pi*f(i,iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)*fv)  )     ) )**gammae2*f(i,iii) 
                            where (g1<0d0)
                                g1=0d0
                            end where
                            g2=LAMBDA1*( 1d0-cc/(cc+ss*(   max(0d0,pi*f(max(i-1,1),iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)*fv)  )       ) )**gammae2*f(max(i-1,1),iii) 
                            where (g2<0d0)
                                g2=0d0
                            end where
                            temp1=pi*f(i,iii)-x(j)*x_vp-g1
                            temp2=pi*f(max(i-1,1),iii)-x(j)*x_vp-g2
                            temp1b=(-cc+(1d0-ss)*(temp1+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)))
                            temp2b=(-cc+(1d0-ss)*(temp2+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)))

                            where (temp1b<1e-10)
                                temp1b=1e-10
                            end where
                            where (temp2b<1e-10)
                                temp2b=1e-10
                            end where

                            vv_l=gacp_l(i,j,ii,jj,:,:)*(kesi*temp1b(:,:)**(1d0-omega)/(1d0-omega)+(1d0-kesi)*temp2b(:,:)**(1d0-omega)/(1d0-omega))&
                                +(1d0-gacp_l(i,j,ii,jj,:,:))*((pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))
                            !vv_record_big_3(i,j,ii,jj,:,:,1)=vv_h
                            !vv_record_big_3(i,j,ii,jj,:,:,2)=vv_m
                            !vv_record_big_3(i,j,ii,jj,:,:,3)=vv_l

                            do inn=1,numi
                                p2(:,:)=p3(:,:,inn)
                                !print *, p2(1,ii),p2(2,ii),p2(3,ii)
                                !pause
                                vv=p2(1,ii)*vv_h+p2(2,ii)*vv_m+p2(3,ii)*vv_l

                                vv_record_big(i,j,ii,jj,:,:,inn)=vv

                                !valone_record(i,j,ii,jj)=(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega)



                                vvc(1,:) =maxval(vv(:,:),1)
                                idxc(1,:)=maxloc(vv(:,:),1)
                                vvs=maxval(vvc(1,:),1)
                                idxs=maxloc(vvc(1,:),1)

                                indexs(i,j,ii,jj)=idxs
                                indexc(i,j,ii,jj)=idxc(1,idxs)
                                !print *,  'indexs(i,j,ii,jj)',indexs(i,j,ii,jj)


                                if (vvs>(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega)) then !(pi*za(i))^(1-omega)/(1-omega)
                                    offermade(i,j,ii,jj)=1
                                else
                                    offermade(i,j,ii,jj)=0
                                end if
                                ! vv_record(i,j,ii,jj,:,:)=vvs

                                Vanew(i,j,ii,jj)=max(vvs,(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))
                                ! vv_record2(i,j,ii,jj,:,:)=vv

                                gc(i,j,ii,jj)=cc(indexc(i,j,ii,jj),indexs(i,j,ii,jj))
                                gs(i,j,ii,jj)=ss(indexc(i,j,ii,jj),indexs(i,j,ii,jj))
                                if (offermade(i,j,ii,jj)==0) then
                                    gc(i,j,ii,jj)=0d0
                                    gs(i,j,ii,jj)=0d0
                                    indexc(i,j,ii,jj)=1
                                    indexs(i,j,ii,jj)=1
                                    vv_record_big(i,j,ii,jj,:,:,inn)=(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega)
                                end if


                                indexs_i(i,j,ii,jj,inn)= indexs(i,j,ii,jj)
                                indexc_i(i,j,ii,jj,inn)= indexc(i,j,ii,jj)
                                gs_i(i,j,ii,jj,inn)= gs(i,j,ii,jj)
                                gc_i(i,j,ii,jj,inn)= gc(i,j,ii,jj)
                                offermade_i(i,j,ii,jj,inn)=offermade(i,j,ii,jj)
                                Vanew_i(i,j,ii,jj,inn)=Vanew(i,j,ii,jj)
                            enddo

                            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                            ! Perfect Information Case
                            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


                            vvs=maxval(vv_m(2,:),1)
                            idxs=maxloc(vv_m(2,:),1)
                            indexs_perfect(i,j,ii,jj)=idxs
                            indexc_perfect(i,j,ii,jj)=2

                            if (vvs>(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega)) then !(pi*za(i))^(1-omega)/(1-omega)
                                offermade_perfect(i,j,ii,jj)=1
                            else
                                offermade_perfect(i,j,ii,jj)=0
                            end if

                            Vanew_perfect(i,j,ii,jj)=max(vvs,(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))

                            gc_perfect(i,j,ii,jj)=cc(indexc_perfect(i,j,ii,jj),indexs_perfect(i,j,ii,jj))
                            gs_perfect(i,j,ii,jj)=ss(indexc_perfect(i,j,ii,jj),indexs_perfect(i,j,ii,jj))
                            if (offermade_perfect(i,j,ii,jj)==0) then
                                gc_perfect(i,j,ii,jj)=0d0
                                gs_perfect(i,j,ii,jj)=0d0
                                indexs_perfect(i,j,ii,jj)=1
                                indexc_perfect(i,j,ii,jj)=1
                            end if




                        end do
                    end do
                end do
            end do

            !$omp end  do
            !$omp end parallel
            !call system_clock(int2)

            !print *, 'Third Part',(int2-int1)

            !call system_clock(int1)






            !now we go back to record target value function value

            !$OMP PARALLEL PRIVATE(i,j,ii,jj,inn,iii,cc,ss,c_value,s_value,temp2s,temp2new,g1s,g2s,temp1h,temp1l,temp1s,temp1_h,temp1_m,temp1_l,x_vp)
            !$OMP DO  SCHEDULE(dynamic)
            do i=1,numz  !acquirer
                do j=1,numx !acquirer
                    do ii=1,numz !target
                        do jj=1,numx

                            ss=reshape(s1,(/3,nums/),s1,(/2,1/))

                            !print *,'ss(1,:)',ss(1,:)
                            !  pause
                            cc=0d0 !zeros(3,nums)

                            cc(1,:) = c_threshold2(i,j,min(ii+1,numz),jj,:)
                            cc(2,:) = c_threshold2(i,j,ii,jj,:)
                            cc(3,:) = c_threshold2(i,j,max(ii-1,1),jj,:)
                            where (cc<0d0)
                                cc=0d0
                            end where

                            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                            ! Imperfect Information Case
                            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                            do inn=1,numi


                                iii=min(ii+1,numz)
                                c_value=gc_i(i,j,iii,jj,inn)
                                s_value=gs_i(i,j,iii,jj,inn)


                                !three cases for target


                                temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                temp2s=max(temp2s,1e-10)  
                                temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                g1s=max(g1s,0d0)
                                g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii) 
                                g2s=max(g2s,0d0)
                                temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j))
                                temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j))
                                temp1h = max(temp1h,1e-10)
                                temp1l = max(temp1l,1e-10)
                                temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                if (offermade_i(i,j,iii,jj,inn)==1) then

                                    temp1_h =  gacp_h(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn))*temp1s+(1d0- gacp_h(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn)))*temp2new
                                else
                                    temp1_h=temp2new
                                end if

                                if (temp1_h > 100000d0) then
                                    print *, 'c_value, s_value, g1s, g2s ', c_value, s_value, g1s, g2s
                                    print *, 'i, j, iii, jj, inn, offermade_i(i,j,iii,jj,inn) ', i, j, iii, jj, inn, offermade_i(i,j,iii,jj,inn)
                                    pause
                                end if


                                iii=ii
                                c_value=gc_i(i,j,iii,jj,inn)
                                s_value=gs_i(i,j,iii,jj,inn)


                                temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                temp2s=max(temp2s,1e-10)
                                temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                g1s=max(g1s,0d0)
                                g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii) 
                                g2s=max(g2s,0d0)
                                temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)) 
                                temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)) 
                                temp1h = max(temp1h,1e-10)
                                temp1l = max(temp1l,1e-10)
                                temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                if (offermade_i(i,j,iii,jj,inn)==1) then

                                    temp1_m =  gacp_m(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn))*temp1s+(1d0- gacp_m(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn)))*temp2new
                                else
                                    temp1_m=temp2new
                                end if

                                if (temp1_m > 100000d0) then
                                    print *, 'c_value, s_value, g1s, g2s ', c_value, s_value, g1s, g2s
                                    print *, 'i, j, iii, jj, inn, offermade_i(i,j,iii,jj,inn) ', i, j, iii, jj, inn, offermade_i(i,j,iii,jj,inn)
                                    pause
                                end if


                                iii=max(ii-1,1)
                                
                                c_value=gc_i(i,j,iii,jj,inn)
                                s_value=gs_i(i,j,iii,jj,inn)



                                temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                temp2s=max(temp2s,1e-10)
                                temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                g1s=max(g1s,0d0)
                                g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii) 
                                g2s=max(g2s,0d0)


                                temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)) 
                                temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)) 
                                temp1h = max(temp1h,1e-10)
                                temp1l = max(temp1l,1e-10)
                                temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                if (offermade_i(i,j,iii,jj,inn)==1) then

                                    temp1_l =  gacp_l(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn))*temp1s+(1d0- gacp_l(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn)))*temp2new
                                else
                                    temp1_l=temp2new
                                end if

                                if (temp1_l > 100000d0) then
                                    print *, 'c_value, s_value, g1s, g2s ', c_value, s_value, g1s, g2s
                                    print *, 'i, j, iii, jj, inn, offermade_i(i,j,iii,jj,inn) ', i, j, iii, jj, inn, offermade_i(i,j,iii,jj,inn)
                                    pause
                                end if


                                Vtnew_h(i,j,ii,jj)=temp1_h
                                Vtnew_m(i,j,ii,jj)=temp1_m
                                Vtnew_l(i,j,ii,jj)=temp1_l


                                Vtnew_h_i(i,j,ii,jj,inn)=temp1_h
                                Vtnew_m_i(i,j,ii,jj,inn)=temp1_m
                                Vtnew_l_i(i,j,ii,jj,inn)=temp1_l

                            enddo


                            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                            ! Perfect Information Case
                            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                            iii=ii
                            c_value=cc(indexc_perfect(i,j,iii,jj),indexs_perfect(i,j,iii,jj))
                            s_value=ss(indexc_perfect(i,j,iii,jj),indexs_perfect(i,j,iii,jj))


                            temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                            temp2s=max(temp2s,1e-10)
                            temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                            x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                            g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                            g1s=max(g1s,0d0)
                            g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii) 
                            g2s=max(g2s,0d0)
                            temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)) 
                            temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)) 
                            temp1h = max(temp1h,1e-10)
                            temp1l = max(temp1l,1e-10)
                            temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                            if (offermade_perfect(i,j,iii,jj)==1) then

                                temp1_m =  gacp_m(i,j,iii,jj,indexc_perfect(i,j,iii,jj),indexs_perfect(i,j,iii,jj))*temp1s+(1d0- gacp_m(i,j,iii,jj,indexc_perfect(i,j,iii,jj),indexs_perfect(i,j,iii,jj)))*temp2new
                            else
                                temp1_m=temp2new
                            end if

                            Vtnew_perfect(i,j,ii,jj) = temp1_m

                        end do
                    end do
                end do
            end do
            !
            !            !
            !$omp end  do
            !$omp end parallel


            !call system_clock(int2)

            !print *, 'Fourth Part',(int2-int1)

            !call system_clock(int1)



            do i=1,numz  ! my pre-innovation productivity
                do j=1,numx ! my merger cost

                    do inn=1,numi ! my innovation choice

                        do ii=1,numz ! other pre-innovation productivity
                            do jj=1,numx ! other merger cost

                                tempt = 0d0
                                do k = -1,1 ! my innovation realization
                                    do m = -1,1 ! other innovation realization
                                        ip = max(min(i+k,numz),1)
                                        iip = max(min(ii+m,numz),1)

                                        if (k == -1) then
                                            tempt = tempt + (info_asym*Vtnew_h_i(iip,jj,ip,j,inn) + (1d0-info_asym)*Vtnew_perfect(iip,jj,ip,j))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                        elseif (k == 0) then
                                            tempt = tempt + (info_asym*Vtnew_m_i(iip,jj,ip,j,inn) + (1d0-info_asym)*Vtnew_perfect(iip,jj,ip,j))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                        else
                                            tempt = tempt + (info_asym*Vtnew_l_i(iip,jj,ip,j,inn) + (1d0-info_asym)*Vtnew_perfect(iip,jj,ip,j))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                        endif
                                    enddo
                                enddo

                                tempa = 0d0
                                do k = -1,1 ! my innovation realization
                                    do m = -1,1 ! other innovation realization
                                        ip = max(min(i+k,numz),1)
                                        iip = max(min(ii+m,numz),1)
                                        tempa = tempa + (info_asym*Vanew_i(ip,j,ii,jj,i_policy(ii,jj)) + (1d0-info_asym)*Vanew_perfect(ip,j,iip,jj))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                    enddo
                                enddo

                                if (i < ii) then
                                    V_meet_i(ii,jj,inn) = tempt
                                elseif (i == ii) then
                                    V_meet_i(ii,jj,inn) = (tempt + tempa)/2d0
                                else
                                    V_meet_i(ii,jj,inn) = tempa
                                endif
                            enddo
                        enddo

                        temp_alone = 0d0
                        do k = -1,1 ! my innovation realization
                            ip = max(min(i+k,numz),1)
                            temp_alone = temp_alone + (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))*p3(2-k,i,inn)
                        enddo

                        Vnew_i(i,j,inn) = eta*(sum(V_meet_i(:,:,inn)*Fs(:,:)/sum(Fs))) + (1d0-eta)*temp_alone - (i_cost(inn)*za(i))

                    enddo


                    if (GE_fixed_Fs_innov==1) then
                        Vnew(i,j) =Vnew_i(i,j,i_policy(i,j))
                    else
                        Vnew(i,j) =maxval(Vnew_i(i,j,:),1)
                        i_policy_new(i,j) =maxloc(Vnew_i(i,j,:),1)
                    end if

                enddo
            enddo

            do i=1,numz
                do j=1,numx
                    v_prod_new(i,j) = 0d0
                    do k = -1,1 ! my innovation realization
                        ip = max(min(i+k,numz),1)
                        v_prod_new(i,j) = v_prod_new(i,j) + (pi*za(ip)+(1d0-psi)/(1d0+r)*v_prod(ip,j))*p2(2-k,i)     
                    enddo
                enddo
            enddo




            diffv=maxval(abs(Vnew-V))
            diff_pol=maxval(abs(i_policy_new - i_policy))

            counter=counter+1



            Va=Vanew
            Vt=Vtnew
            V=Vnew*vfi_update+ V*(1d0-vfi_update)
            v_prod = v_prod_new
            i_policy = i_policy_new





        end do




        OPEN(1,file='V_record.txt') ; WRITE(1,"(f20.8)") V_record


        iter30=0
        diff_dist_inner=1d0
        if (diff_dist<1e-8) then

            maxiter30=1
        else
            maxiter30=maxiter30_value
        end if




        do while (iter30<maxiter30.and.diff_dist_inner>1e-8)




            gacp_eqm = 0
            gacp_eqm_perfect = 0
            do i = 1,numz
                do j = 1,numx
                    do ii = 1,numz

                        do jj = 1,numx
                            do inn = 1,numi
                                gacp_eqm_i(i,j,ii,jj,1,inn) = nint( gacp_h(i,j,ii,jj,indexc_i(i,j,ii,jj,inn),indexs_i(i,j,ii,jj,inn)))
                                gacp_eqm_i(i,j,ii,jj,2,inn) = nint( gacp_m(i,j,ii,jj,indexc_i(i,j,ii,jj,inn),indexs_i(i,j,ii,jj,inn)))
                                gacp_eqm_i(i,j,ii,jj,3,inn) = nint( gacp_l(i,j,ii,jj,indexc_i(i,j,ii,jj,inn),indexs_i(i,j,ii,jj,inn)))
                            enddo
                            gacp_eqm_perfect(i,j,ii,jj) = nint( gacp_m(i,j,ii,jj,indexc_perfect(i,j,ii,jj),indexs_perfect(i,j,ii,jj)))
                        end do
                    end do
                end do
            end do



            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NEW TRANS CODE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            trans = 0d0
            dummy = 0d0

            do j = 1,numx ! x
                do i = 1,numz ! z
                    ! exogenous exit
                    do ii_e = 1,numz ! entrant z
                        do jj_e = 1,numx ! entrant x
                            call giveindex(i,j)
                            a=zx_idx !today index
                            call giveindex(ii_e,jj_e)
                            b=zx_idx !tomorrow index

                            F_exog_corrected = F_exog(ii_e,jj_e)/sum(F_exog)
                            if (sum(F_exog)==0d0) then
                                F_exog_corrected = 0d0
                            end if

                            trans(a,b) = trans(a,b) + psi*F_exog_corrected
                            dummy(i,j) = dummy(i,j) +  psi*F_exog_corrected
                        end do
                    end do
                    ! no exit
                    ! no meeting
                    ! high z'
                    call giveindex(i,j)
                    a=zx_idx !today index
                    call giveindex(min(i+1,numz),j)
                    b=zx_idx !tomorrow index
                    trans(a,b) = trans(a,b) + (1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
                    dummy(i,j) = dummy(i,j) + (1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
                    ! med z'
                    call giveindex(i,j)
                    a=zx_idx !today index
                    call giveindex(i,j)
                    b=zx_idx !tomorrow index
                    trans(a,b) = trans(a,b) + (1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
                    dummy(i,j) = dummy(i,j) + (1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
                    ! low z'
                    call giveindex(i,j)
                    a=zx_idx !today index
                    call giveindex(max(i-1,1),j)
                    b=zx_idx !tomorrow index
                    trans(a,b) = trans(a,b) + (1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)
                    dummy(i,j) = dummy(i,j) + (1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)


                    !if (be_acq(i,j) == 1) then ! Acquirer
                    do ii = 1,numz ! other pre-innovation productivity
                        do jj = 1,numx ! other fixed cost

                            if (i>ii) then
                                prob_acq(i,j) = 1d0
                            elseif (i==ii) then
                                prob_acq(i,j) = 0.5
                            else
                                prob_acq(i,j) = 0d0
                            endif

                            ! meet, acquirer -- imperfect information case
                            do k = -1,1 ! acquirer innovation realization
                                do m = -1,1 ! target innovation realization

                                    ip = max(min(i+k,numz),1)
                                    iip = max(min(ii+m,numz),1)
                                    tttr = Fs(ii,jj)/sum(Fs)*info_asym*prob_acq(i,j) *eta*p3(2-k,i,i_policy(i,j))*p3(2-m,ii,i_policy(ii,jj))*(1d0-psi)
                                    if (sum(Fs)==0d0) then
                                        tttr = 0d0
                                    end if
                                    if (offermade_i(ip,j,ii,jj,i_policy(ii,jj)) == 1) then


                                        if (gacp_eqm_i(ip,j,ii,jj,2-m,i_policy(ii,jj)) == 1) then! if target accepts
                                            ! merger risk -- no drop
                                            call zinterp(ip,iip,p_lb, p_ub,i_lb,i_ub)


                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(i_lb,j)
                                            b=zx_idx !tomorrow index



                                            trans(a,b)= trans(a,b)+ tttr*kesi*p_lb
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(i_ub,j)
                                            b=zx_idx !tomorrow index


                                            trans(a,b)= trans(a,b)+ kesi*p_ub*tttr
                                            dummy(i,j) = dummy(i,j) + kesi*p_lb*tttr
                                            dummy(i,j) = dummy(i,j) + kesi*p_ub*tttr

                                            ! merger risk -- drop
                                            !   [ p_lb, p_ub,i_lb,i_ub]=zinterp(max(ip-1,1),iip)

                                            call zinterp(max(ip-1,1),iip,p_lb, p_ub,i_lb,i_ub)
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(i_lb,j)
                                            b=zx_idx !tomorrow index

                                            trans(a,b)= trans(a,b)+ (1d0-kesi)*p_lb*tttr
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(i_ub,j)
                                            b=zx_idx !tomorrow index


                                            trans(a,b)= trans(a,b)+ (1d0-kesi)*p_ub*tttr
                                            dummy(i,j) = dummy(i,j) +  (1d0-kesi)*p_lb*tttr
                                            dummy(i,j) = dummy(i,j) +(1d0-kesi)*p_ub*tttr

                                        else ! if target rejects


                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call  giveindex(ip,j)
                                            b=zx_idx !tomorrow index


                                            trans(a,b) = trans(a,b) + tttr
                                            dummy(i,j) = dummy(i,j) +  tttr
                                        end if
                                    else ! no offer
                                        call giveindex(i,j)
                                        a=zx_idx !today index
                                        call  giveindex(ip,j)
                                        b=zx_idx !tomorrow index


                                        trans(a,b) = trans(a,b) + tttr
                                        dummy(i,j) = dummy(i,j) + tttr
                                    end if
                                end do
                            end do
                        end do
                    end do

                    ! meet -- perfect information case
                    do ii = 1,numz ! other pre-innovation productivity
                        do jj = 1,numx ! other fixed cost

                            if (i>ii) then
                                prob_acq(i,j) = 1d0
                            elseif (i==ii) then
                                prob_acq(i,j) = 0.5
                            else
                                prob_acq(i,j) = 0d0
                            endif

                            do k = -1,1 ! acquirer innovation realization
                                do m = -1,1 ! target innovation realization
                                    ip = max(min(i+k,numz),1)
                                    iip = max(min(ii+m,numz),1)
                                    ! tttr = Fsp(iip,jj)/sum(Fsp)*(1d0-info_asym)*prob_acq(i,j)
                                    ! if (sum(Fsp)==0d0) then
                                    !     tttr = 0d0
                                    ! end if
                                    tttr = Fs(ii,jj)/sum(Fs)*(1d0-info_asym)*prob_acq(i,j) *eta*p3(2-k,i,i_policy(i,j))*p3(2-m,ii,i_policy(ii,jj))*(1d0-psi)
                                    if (sum(Fs)==0d0) then
                                        tttr = 0d0
                                    end if
                                    if (offermade_perfect(ip,j,iip,jj) == 1) then


                                        if (gacp_eqm_perfect(ip,j,iip,jj) == 1) then! if target accepts
                                            ! merger risk -- no drop
                                            call zinterp(ip,iip,p_lb, p_ub,i_lb,i_ub)



                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(i_lb,j)
                                            b=zx_idx !tomorrow index



                                            trans(a,b)= trans(a,b)+ kesi*p_lb*tttr
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(i_ub,j)
                                            b=zx_idx !tomorrow index


                                            trans(a,b)= trans(a,b)+ kesi*p_ub*tttr
                                            dummy(i,j) = dummy(i,j) + kesi*p_lb*tttr
                                            dummy(i,j) = dummy(i,j) + kesi*p_ub*tttr

                                            ! merger risk -- drop
                                            !   [ p_lb, p_ub,i_lb,i_ub]=zinterp(max(ip-1,1),iip)

                                            call zinterp(max(ip-1,1),iip,p_lb, p_ub,i_lb,i_ub)
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(i_lb,j)
                                            b=zx_idx !tomorrow index

                                            trans(a,b)= trans(a,b)+ (1d0-kesi)*p_lb*tttr
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(i_ub,j)
                                            b=zx_idx !tomorrow index



                                            trans(a,b)= trans(a,b)+ (1d0-kesi)*p_ub*tttr
                                            dummy(i,j) = dummy(i,j) +  (1d0-kesi)*p_lb*tttr
                                            dummy(i,j) = dummy(i,j) + (1d0-kesi)*p_ub*tttr

                                        else ! if target rejects


                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call  giveindex(ip,j)
                                            b=zx_idx !tomorrow index


                                            trans(a,b) = trans(a,b) + tttr
                                            dummy(i,j) = dummy(i,j) +  tttr
                                        end if
                                    else ! no offer
                                        call giveindex(i,j)
                                        a=zx_idx !today index
                                        call  giveindex(ip,j)
                                        b=zx_idx !tomorrow index



                                        trans(a,b) = trans(a,b) + tttr
                                        dummy(i,j) = dummy(i,j) + tttr
                                    end if
                                end do
                            end do
                        end do
                    end do

                    ! meet -- imperfect information case
                    do ii = 1,numz ! other pre-innovaton productivity
                        do jj = 1,numx ! other fixed cost

                            if (i>ii) then
                                prob_acq(i,j) = 1d0
                            elseif (i==ii) then
                                prob_acq(i,j) = 0.5
                            else
                                prob_acq(i,j) = 0d0
                            endif

                            do k = -1,1 ! my innovation realization
                                do m = -1,1 ! other innovation realization

                                    ip = max(min(i+k,numz),1) ! target after innov
                                    iip = max(min(ii+m,numz),1) ! acquirer after innov

                                    aaa = Fs(ii,jj)/sum(Fs)*info_asym*(1d0-prob_acq(i,j)) *eta*p3(2-k,i,i_policy(i,j))*p3(2-m,ii,i_policy(ii,jj))*(1d0-psi)
                                    if (sum(Fs) == 0d0) then
                                        aaa = 0d0
                                    end if
                                    if (offermade_i(iip,jj,i,j,i_policy(i,j)) == 1) then
                                        if (gacp_eqm_i(iip,jj,i,j,2-k,i_policy(i,j)) == 1) then ! if target accepts
                                            do ii_e = 1,numz ! entrant z
                                                do jj_e = 1,numx ! entrant x
                                                    call giveindex(i,j)
                                                    a=zx_idx !today index
                                                    call giveindex(ii_e,jj_e)
                                                    b=zx_idx

                                                    F_exog_corrected = F_exog(ii_e,jj_e)/sum(F_exog)
                                                    if (sum(F_exog)==0d0) then
                                                        F_exog_corrected = 0d0
                                                    end if


                                                    trans(a,b) = trans(a,b) + aaa*(F_exog_corrected)
                                                    dummy(i,j) = dummy(i,j) +  aaa*(F_exog_corrected)
                                                end do
                                            end do
                                        else ! if target rejects
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(ip,j)
                                            b=zx_idx



                                            trans(a,b) = trans(a,b) + aaa
                                            dummy(i,j) = dummy(i,j)  +aaa
                                        end if
                                    else ! no offer made
                                        call giveindex(i,j)
                                        a=zx_idx !today index
                                        call giveindex(ip,j)
                                        b=zx_idx


                                        trans(a,b) = trans(a,b) + aaa
                                        dummy(i,j) = dummy(i,j)  + aaa
                                    end if
                                end do
                            end do
                        end do
                    end do


                    ! meet -- perfect information case
                    do ii = 1,numz ! other pre-innovaton productivity
                        do jj = 1,numx ! other fixed cost

                            if (i>ii) then
                                prob_acq(i,j) = 1d0
                            elseif (i==ii) then
                                prob_acq(i,j) = 0.5
                            else
                                prob_acq(i,j) = 0d0
                            endif

                            do k = -1,1 ! my innovation realization
                                do m = -1,1 ! other innovation realization

                                    ip = max(min(i+k,numz),1) ! target after innov
                                    iip = max(min(ii+m,numz),1) ! acquirer after innov

                                    aaa = Fs(ii,jj)/sum(Fs)*(1d0-info_asym)*(1d0-prob_acq(i,j)) *eta*p3(2-k,i,i_policy(i,j))*p3(2-m,ii,i_policy(ii,jj))*(1d0-psi)
                                    if (sum(Fs) == 0d0) then
                                        aaa = 0d0
                                    end if

                                    if (offermade_perfect(iip,jj,ip,j) == 1) then
                                        if (gacp_eqm_perfect(iip,jj,ip,j) == 1) then ! if target accepts
                                            do ii_e = 1,numz ! entrant z
                                                do jj_e = 1,numx ! entrant x
                                                    call giveindex(i,j)
                                                    a=zx_idx !today index
                                                    call giveindex(ii_e,jj_e)
                                                    b=zx_idx

                                                    F_exog_corrected = F_exog(ii_e,jj_e)/sum(F_exog)
                                                    if (sum(F_exog)==0d0) then
                                                        F_exog_corrected = 0d0
                                                    end if

                                                    trans(a,b) = trans(a,b) + aaa*(F_exog_corrected)
                                                    dummy(i,j) = dummy(i,j) + aaa*(F_exog_corrected)
                                                end do
                                            end do
                                        else ! if target rejects
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(ip,j)
                                            b=zx_idx

                                            trans(a,b) = trans(a,b) + aaa
                                            dummy(i,j) = dummy(i,j)  + aaa
                                        end if
                                    else ! no offer made
                                        call giveindex(i,j)
                                        a=zx_idx !today index
                                        call giveindex(ip,j)
                                        b=zx_idx
                                        trans(a,b) = trans(a,b) + aaa
                                        dummy(i,j) = dummy(i,j)  + aaa
                                    end if
                                end do
                            end do
                        end do
                    end do
                    !end if
                end do
            end do

            stat_dist=0d0

            do i = 1,numz
                do j = 1,numx
                    call giveindex(i,j)
                    b=zx_idx

                    stat_dist(b,1) = Fs(i,j)
                end do
            end do


            stat_dist_new = matmul(transpose(trans),stat_dist)


            Fs_old = Fs
            if (GE_fixed_Fs_innov==1) then

                Fs=Fs_load
                Fs_old=Fs_load
            else

                do i = 1,numz
                    do j = 1,numx
                        call giveindex(i,j)
                        b=zx_idx
                        Fs(i,j) = stat_dist_new(b,1)
                    end do
                end do
            end if

            diff_dist_inner=maxval(abs(Fs_old-Fs))

            iter30=iter30+1


        end do


        iter_dist = iter_dist + 1


        if (iter_dist==1) then
            diff_dist = 1000d0
        else
            diff_dist=sqrt(sum((Fs_record-Fs)**2d0))/(1d0+sqrt((sum(Fs_record**2d0))))

        endif


        Fs_record = Fs

        zbar = 0d0
        do i=1,numz
            do j=1,numx
                zbar = zbar + zt(i)*Fs(i,j)
            enddo
        enddo


        print *,'iter_dist diff_dist (distribution)', iter_dist, diff_dist





    end do




    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if (illustration==1) then




        !acquirer problem

        allocate(vv_record_big_f(numz,numx,numz,numx,numc,nums,numi))
        allocate(vv_record_big_f_short(numz,numx,numz,numx,numc,nums))
        allocate( gacp_h_f(numz,numx,numz,numx,numc,nums))
        allocate(gacp_m_f(numz,numx,numz,numx,numc,nums))
        allocate(gacp_l_f(numz,numx,numz,numx,numc,nums))
        allocate(Vtnew_h_f(numz,numx,numz,numx,numc,nums,numi))
        allocate(Vtnew_m_f(numz,numx,numz,numx,numc,nums,numi))
        allocate(Vtnew_l_f(numz,numx,numz,numx,numc,nums,numi))
        allocate(Vtnew_h_f_short(numz,numx,numz,numx,numc,nums))
        allocate(Vtnew_m_f_short(numz,numx,numz,numx,numc,nums))
        allocate(Vtnew_l_f_short(numz,numx,numz,numx,numc,nums))

        !print *, ss_f(1,1:nums)
        !print *, s1
        !pause


        vv_record_big_f_short=0d0
        vv_record_big_f=0d0
        gacp_h_f=0d0
        gacp_m_f=0d0
        gacp_l_f=0d0
        c1(1,16)=30.470598130000000  

        ss_f=reshape(s1,(/numc,nums/),s1,(/2,1/))
        cc_f=reshape(c1,(/numc,nums/),c1,(/1,2/))


        !$OMP PARALLEL PRIVATE(i,j,ii,jj,inn,iii,cc,cc1_f,cc2_f,cc3_f,g1_f,g2_f,temp1_f,temp2_f,temp1b_f,temp2b_f,vv_h_f,vv_m_f,vv_l_f,vv_f,vvc_f,vvs_f,idxc_f,idxs_f,p2,x_vp)
        !$OMP DO collapse(4) SCHEDULE(dynamic)
        do i=1,numz  !acquirer
            do j=1,numx !acquirer
                do ii=1,numz !observed zt_hat
                    do jj=1,numx


                        cc=0d0

                        cc(1,:) = c_threshold2(i,j,min(ii+1,numz),jj,:)
                        cc(2,:) = c_threshold2(i,j,ii,jj,:)
                        cc(3,:) = c_threshold2(i,j,max(ii-1,1),jj,:)
                        where (cc<0d0)
                            cc=0d0
                        end where


                        cc1_f=reshape(cc(1,:),(/numc,nums/),cc(1,:),(/2,1/))
                        cc2_f=reshape(cc(2,:),(/numc,nums/),cc(2,:),(/2,1/))
                        cc3_f=reshape(cc(3,:),(/numc,nums/),cc(3,:),(/2,1/))

                        where (cc_f>=cc1_f)
                            gacp_h_f(i,j,ii,jj,:,:)=1d0  
                        elsewhere
                            gacp_h_f(i,j,ii,jj,:,:)=0d0  
                        end where

                        where (cc_f>=cc2_f)
                            gacp_m_f(i,j,ii,jj,:,:)=1d0
                        elsewhere
                            gacp_m_f(i,j,ii,jj,:,:)=0d0  

                        end where

                        where (cc_f>=cc3_f)
                            gacp_l_f(i,j,ii,jj,:,:)=1d0
                        elsewhere
                            gacp_l_f(i,j,ii,jj,:,:)=0d0  

                        end where

                        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                        ! Imperfect Information Case
                        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                        iii = min(ii+1,numz)
                        x_vp = ((zt(iii)/zt(i))**phi * (zt(iii)+zt(i))**phi2)
                        g1_f=LAMBDA1*( 1d0-cc_f/(cc_f+ss_f*(  max(0d0,pi*f(i,iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)*fv)  )     ) )**gammae2*f(i,iii) 


                        where (g1_f<0d0)
                            g1_f=0d0
                        end where




                        g2_f=LAMBDA1*( 1d0-cc_f/(cc_f+ss_f*(   max(0d0,pi*f(max(i-1,1),iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)*fv)  )       ) )**gammae2*f(max(i-1,1),iii) 
                        where (g2_f<0d0)
                            g2_f=0d0
                        end where

                        temp1_f=pi*f(i,iii)-x(j)*x_vp-g1_f
                        temp2_f=pi*f(max(i-1,1),iii)-x(j)*x_vp-g2_f
                        temp1b_f=(-cc_f+(1d0-ss_f)*(temp1_f+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)))
                        temp2b_f=(-cc_f+(1d0-ss_f)*(temp2_f+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)))


                        where (temp1b_f<1e-10)
                            temp1b_f=1e-10
                        end where
                        where (temp2b_f<1e-10)
                            temp2b_f=1e-10
                        end where


                        vv_h_f=gacp_h_f(i,j,ii,jj,:,:)*(kesi*temp1b_f(:,:)**(1d0-omega)/(1d0-omega)+(1d0-kesi)*temp2b_f(:,:)**(1d0-omega)/(1d0-omega))&
                            +(1d0- gacp_h_f(i,j,ii,jj,:,:))*((pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))



                        iii=ii
                        x_vp = ((zt(iii)/zt(i))**phi * (zt(iii)+zt(i))**phi2)
                        g1_f=LAMBDA1*( 1d0-cc_f/(cc_f+ss_f*(  max(0d0,pi*f(i,iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)*fv)  )     ) )**gammae2*f(i,iii) 
                        where (g1_f<0d0)
                            g1_f=0d0
                        end where
                        g2_f=LAMBDA1*( 1d0-cc_f/(cc_f+ss_f*(   max(0d0,pi*f(max(i-1,1),iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)*fv)  )       ) )**gammae2*f(max(i-1,1),iii) 
                        where (g2_f<0d0)
                            g2_f=0d0
                        end where
                        temp1_f=pi*f(i,iii)-x(j)*x_vp-g1_f
                        temp2_f=pi*f(max(i-1,1),iii)-x(j)*x_vp-g2_f
                        temp1b_f=(-cc_f+(1d0-ss_f)*(temp1_f+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)))
                        temp2b_f=(-cc_f+(1d0-ss_f)*(temp2_f+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)))

                        where (temp1b_f<1e-10)
                            temp1b_f=1e-10
                        end where
                        where (temp2b_f<1e-10)
                            temp2b_f=1e-10
                        end where

                        vv_m_f=gacp_m_f(i,j,ii,jj,:,:)*(kesi*temp1b_f(:,:)**(1d0-omega)/(1d0-omega)+(1d0-kesi)*temp2b_f(:,:)**(1d0-omega)/(1d0-omega))&
                            +(1d0- gacp_m_f(i,j,ii,jj,:,:))*((pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))



                        iii = max(ii-1,1)
                        x_vp = ((zt(iii)/zt(i))**phi * (zt(iii)+zt(i))**phi2)
                        g1_f=LAMBDA1*( 1d0-cc_f/(cc_f+ss_f*(  max(0d0,pi*f(i,iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)*fv)  )     ) )**gammae2*f(i,iii) 
                        where (g1_f<0d0)
                            g1_f=0d0
                        end where
                        g2_f=LAMBDA1*( 1d0-cc_f/(cc_f+ss_f*(   max(0d0,pi*f(max(i-1,1),iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)*fv)  )       ) )**gammae2*f(max(i-1,1),iii) 
                        where (g2_f<0d0)
                            g2_f=0d0
                        end where
                        temp1_f=pi*f(i,iii)-x(j)*x_vp-g1_f
                        temp2_f=pi*f(max(i-1,1),iii)-x(j)*x_vp-g2_f
                        temp1b_f=(-cc_f+(1d0-ss_f)*(temp1_f+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)))
                        temp2b_f=(-cc_f+(1d0-ss_f)*(temp2_f+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)))

                        where (temp1b_f<1e-10)
                            temp1b_f=1e-10
                        end where
                        where (temp2b_f<1e-10)
                            temp2b_f=1e-10
                        end where

                        vv_l_f=gacp_l_f(i,j,ii,jj,:,:)*(kesi*temp1b_f(:,:)**(1d0-omega)/(1d0-omega)+(1d0-kesi)*temp2b_f(:,:)**(1d0-omega)/(1d0-omega))&
                            +(1d0-gacp_l_f(i,j,ii,jj,:,:))*((pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))


                        do inn=1,numi
                            p2(:,:)=p3(:,:,inn)
                            vv_f=p2(1,ii)*vv_h_f+p2(2,ii)*vv_m_f+p2(3,ii)*vv_l_f

                            vv_record_big_f(i,j,ii,jj,:,:,inn)=vv_f



                            vvc_f(1,:) =maxval(vv_f(:,:),1)
                            idxc_f(1,:)=maxloc(vv_f(:,:),1)
                            vvs_f=maxval(vvc_f(1,:),1)
                            idxs_f=maxloc(vvc_f(1,:),1)

                            indexs_f(i,j,ii,jj)=idxs_f
                            indexc_f(i,j,ii,jj)=idxc_f(1,idxs_f)


                            if (vvs_f>(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega)) then !(pi*za(i))^(1-omega)/(1-omega)
                                offermade_f(i,j,ii,jj)=1
                            else
                                offermade_f(i,j,ii,jj)=0
                            end if

                            Vanew_f(i,j,ii,jj)=max(vvs_f,(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))


                            gc_f(i,j,ii,jj)=cc_f(indexc_f(i,j,ii,jj),indexs_f(i,j,ii,jj))
                            gs_f(i,j,ii,jj)=ss_f(indexc_f(i,j,ii,jj),indexs_f(i,j,ii,jj))
                            if (offermade_f(i,j,ii,jj)==0) then
                                gc_f(i,j,ii,jj)=0d0
                                gs_f(i,j,ii,jj)=0d0
                                indexc_f(i,j,ii,jj)=1
                                indexs_f(i,j,ii,jj)=1
                                vv_record_big_f(i,j,ii,jj,:,:,inn)=(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega)  
                            end if


                            indexs_i_f(i,j,ii,jj,inn)= indexs_f(i,j,ii,jj)
                            indexc_i_f(i,j,ii,jj,inn)= indexc_f(i,j,ii,jj)
                            gs_i_f(i,j,ii,jj,inn)= gs_f(i,j,ii,jj)
                            gc_i_f(i,j,ii,jj,inn)= gc_f(i,j,ii,jj)
                            offermade_i_f(i,j,ii,jj,inn)=offermade_f(i,j,ii,jj)
                            Vanew_i_f(i,j,ii,jj,inn)=Vanew_f(i,j,ii,jj)
                        enddo
                        vv_record_big_f_short(i,j,ii,jj,:,:)= vv_record_big_f(i,j,ii,jj,:,:,i_policy(ii,jj))


                    end do
                end do
            end do
        end do

        !$omp end  do
        !$omp end parallel



        !now we go back to record target value function value

        !$OMP PARALLEL PRIVATE(i,j,ii,jj,ic,is,inn,iii,cc,ss,c_value,s_value,temp2s,temp2new,g1s,g2s,temp1h,temp1l,temp1s,temp1_h,temp1_m,temp1_l,x_vp)
        !$OMP DO  SCHEDULE(dynamic)
        do i=1,numz  !acquirer
            do j=1,numx !acquirer
                do ii=1,numz !target
                    do jj=1,numx
                        do ic=1,numc
                            do is=1,nums


                                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                                ! Imperfect Information Case
                                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                                do inn=1,numi


                                    iii=min(ii+1,numz)
                                    c_value=c1(1,ic)
                                    s_value=s1(1,is)

                                    !three cases for target


                                    temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                    temp2s=max(temp2s,1e-10)
                                    temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                    x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                    g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                    g1s=max(g1s,0d0)
                                    g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii) 
                                    g2s=max(g2s,0d0)
                                    temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)) 
                                    temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j))
                                    temp1h = max(temp1h,1e-10)
                                    temp1l = max(temp1l,1e-10)
                                    temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                    if (offermade_i_f(i,j,iii,jj,inn)==1) then

                                        temp1_h =  gacp_h_f(i,j,iii,jj,ic,is)*temp1s+(1d0- gacp_h_f(i,j,iii,jj,ic,is))*temp2new
                                    else
                                        temp1_h=temp2new
                                    end if



                                    iii=ii
                                    c_value=c1(1,ic)
                                    s_value=s1(1,is)



                                    temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                    temp2s=max(temp2s,1e-10)
                                    temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                    x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                    g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                    g1s=max(g1s,0d0)
                                    g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii)
                                    g2s=max(g2s,0d0)
                                    temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j))
                                    temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j))
                                    temp1h = max(temp1h,1e-10)
                                    temp1l = max(temp1l,1e-10)
                                    temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                    if (offermade_i_f(i,j,iii,jj,inn)==1) then

                                        temp1_m =  gacp_m_f(i,j,iii,jj,ic,is)*temp1s+(1d0- gacp_m_f(i,j,iii,jj,ic,is))*temp2new
                                    else
                                        temp1_m=temp2new
                                    end if



                                    iii=max(ii-1,1)
                                    c_value=c1(1,ic)
                                    s_value=s1(1,is)



                                    temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                    temp2s=max(temp2s,1e-10)
                                    temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                    x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                    g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                    g1s=max(g1s,0d0)
                                    g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii)
                                    g2s=max(g2s,0d0)
                                    temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j))
                                    temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j))
                                    temp1h = max(temp1h,1e-10)
                                    temp1l = max(temp1l,1e-10)
                                    temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                    if (offermade_i_f(i,j,iii,jj,inn)==1) then

                                        temp1_l =  gacp_l_f(i,j,iii,jj,ic,is)*temp1s+(1d0- gacp_l_f(i,j,iii,jj,ic,is))*temp2new
                                    else
                                        temp1_l=temp2new
                                    end if



                                    Vtnew_h_f(i,j,ii,jj,ic,is,inn)=temp1_h
                                    Vtnew_m_f(i,j,ii,jj,ic,is,inn)=temp1_m
                                    Vtnew_l_f(i,j,ii,jj,ic,is,inn)=temp1_l



                                enddo


                                Vtnew_h_f_short(i,j,ii,jj,ic,is)=Vtnew_h_f(i,j,ii,jj,ic,is,i_policy(ii,jj))
                                Vtnew_m_f_short(i,j,ii,jj,ic,is)=Vtnew_m_f(i,j,ii,jj,ic,is,i_policy(ii,jj))
                                Vtnew_l_f_short(i,j,ii,jj,ic,is)=Vtnew_l_f(i,j,ii,jj,ic,is,i_policy(ii,jj))
                            end do
                        end do
                    end do
                end do
            end do
        end do
        !
        !            !
        !$omp end  do
        !$omp end parallel









        if (estimate==0 .and. comparastat==0) then
            OPEN(1,file='vv_record_big_f.txt') ; WRITE(1,"(f20.8)") vv_record_big_f
            OPEN(1,file='vv_record_big_f_short.txt') ; WRITE(1,"(f20.8)") vv_record_big_f_short
            OPEN(1,file='c1.txt') ; WRITE(1,"(f20.8)") c1
            OPEN(1,file='s1.txt') ; WRITE(1,"(f20.8)") s1
            OPEN(1,file='gs_i_f.txt') ; WRITE(1,"(f20.8)") gs_i_f
            OPEN(1,file='gc_i_f.txt') ; WRITE(1,"(f20.8)") gc_i_f
            OPEN(1,file='gacp_h_f.txt') ; WRITE(1,"(f20.8)") gacp_h_f
            OPEN(1,file='gacp_m_f.txt') ; WRITE(1,"(f20.8)") gacp_m_f
            OPEN(1,file='gacp_l_f.txt') ; WRITE(1,"(f20.8)") gacp_l_f
            OPEN(1,file='Vtnew_h_f_short.txt') ; WRITE(1,"(f20.8)") Vtnew_h_f_short
            OPEN(1,file='Vtnew_m_f_short.txt') ; WRITE(1,"(f20.8)") Vtnew_m_f_short
            OPEN(1,file='Vtnew_l_f_short.txt') ; WRITE(1,"(f20.8)") Vtnew_l_f_short
        end if
        deallocate(vv_record_big_f)
        deallocate(gacp_h_f)
        deallocate(gacp_m_f)
        deallocate(gacp_l_f)
        deallocate(vv_record_big_f_short)
        deallocate(Vtnew_h_f)
        deallocate(Vtnew_m_f)
        deallocate(Vtnew_l_f)
        deallocate(Vtnew_h_f_short)
        deallocate(Vtnew_m_f_short)
        deallocate(Vtnew_l_f_short)


    end if







    if (estimate==0 .and. comparastat==0) then
        OPEN(1,file='za.txt') ; WRITE(1,"(f20.8)")  za
        OPEN(1,file='zt.txt') ; WRITE(1,"(f20.8)")  zt
        OPEN(1,file='s.txt') ; WRITE(1,"(f20.8)")  s
        OPEN(1,file='f.txt') ; WRITE(1,"(f20.8)")  f

        OPEN(1,file='Fs.txt') ;   WRITE(1,"(f20.8)")  Fs
        OPEN(1,file='stat_dist.txt');    WRITE(1,"(f20.8)")  stat_dist
        OPEN(1,file='trans.txt');    WRITE(1,"(f20.8)")  trans
        OPEN(1,file='gc_i.txt') ;   WRITE(1,"(f20.8)")  gc_i
        OPEN(1,file='gs_i.txt') ;   WRITE(1,"(f20.8)")  gs_i
        OPEN(1,file='indexs_i.txt') ;   WRITE(1,"(I4)")  indexs_i
        OPEN(1,file='indexc_i.txt') ;   WRITE(1,"(I4)")  indexc_i
        OPEN(1,file='i_policy.txt') ;   WRITE(1,"(I4)")  i_policy
        OPEN(1,file='para.txt') ; WRITE(1,"(f20.8)") params
        OPEN(1,file='c_threshold2.txt') ; WRITE(1,"(f20.8)") c_threshold2
        OPEN(1,file='error_c2.txt') ; WRITE(1,"(f20.8)") error_c2

        OPEN(1,file='be_acq.txt') ; WRITE(1,"(I4)") be_acq
        OPEN(1,file='p2.txt') ; WRITE(1,"(f20.8)") p2
        OPEN(1,file='p3.txt') ; WRITE(1,"(f20.8)") p3
        OPEN(1,file='pa.txt') ; WRITE(1,*) pa
        OPEN(1,file='pt.txt') ; WRITE(1,*) pt
        OPEN(1,file='theta.txt') ; WRITE(1,*) theta

        OPEN(1,file='Vtnew_h_i.txt') ; WRITE(1,"(f20.8)") Vtnew_h_i
        OPEN(1,file='Vtnew_m_i.txt') ; WRITE(1,"(f20.8)") Vtnew_m_i
        OPEN(1,file='Vtnew_l_i.txt') ; WRITE(1,"(f20.8)") Vtnew_l_i
        OPEN(1,file='Vanew_i.txt') ; WRITE(1,"(f20.8)") Vanew_i
        OPEN(1,file='V.txt') ; WRITE(1,"(f20.8)") V
        OPEN(1,file='offermade_i.txt') ; WRITE(1,"(I4)") offermade_i
        OPEN(1,file='gacp_h.txt') ; WRITE(1,"(f20.8)") gacp_h
        OPEN(1,file='gacp_m.txt') ; WRITE(1,"(f20.8)") gacp_m
        OPEN(1,file='gacp_l.txt') ; WRITE(1,"(f20.8)") gacp_l
        OPEN(1,file='gacp_eqm_i.txt') ; WRITE(1,"(I4)") gacp_eqm_i
        OPEN(1,file='pi.txt') ; WRITE(1,"(f20.8)") pi
        OPEN(1,file='vv_record2.txt') ; WRITE(1,"(f20.8)") vv_record2
        !OPEN(1,file='vv_record.txt') ; WRITE(1,"(f20.8)") vv_record
        OPEN(1,file='vv_record_big_3.txt') ; WRITE(1,"(f20.8)") vv_record_big_3
        OPEN(1,file='vv_record_big.txt') ; WRITE(1,"(f20.8)") vv_record_big
        OPEN(1,file='valone_record.txt') ; WRITE(1,"(f20.8)") valone_record
        OPEN(1,file='V_interp.txt') ; WRITE(1,"(f20.8)") V_interp
        OPEN(1,file='ss.txt') ; WRITE(1,"(f20.8)") ss

        OPEN(1,file='gc_perfect.txt') ; WRITE(1,"(f20.8)") gc_perfect
        OPEN(1,file='gs_perfect.txt') ; WRITE(1,"(f20.8)") gs_perfect
        OPEN(1,file='Vanew_perfect.txt') ; WRITE(1,"(f20.8)") Vanew_perfect
        OPEN(1,file='Vtnew_perfect.txt') ; WRITE(1,"(f20.8)") Vtnew_perfect
        OPEN(1,file='Ftp.txt') ; WRITE(1,"(f20.8)") Ftp

        OPEN(1,file='indexs_perfect.txt') ;   WRITE(1,"(I4)")  indexs_perfect
        OPEN(1,file='indexc_perfect.txt') ;   WRITE(1,"(I4)")  indexc_perfect
        OPEN(1,file='offermade_perfect.txt') ;   WRITE(1,"(I4)")  offermade_perfect
        OPEN(1,file='gacp_eqm_perfect.txt') ;   WRITE(1,"(I4)")  gacp_eqm_perfect
        OPEN(1,file='prob_acq.txt') ;   WRITE(1,"(f20.8)")  prob_acq
        OPEN(1,file='dist_improve.txt') ;   WRITE(1,"(f20.8)")  dist_improve
        OPEN(1,file='dist_record.txt') ;   WRITE(1,"(f20.8)")  dist_record



        !print *,valone_record(:,:,1,1)

    end if



    deallocate(V_record)






    end subroutine policy



    subroutine makemoments(modelmom2)

    USE parameters
    USE globals
    use tools



    IMPLICIT NONE

    include "mkl.fi"

    real*8::corr_equity_zt,mean_equity_ratio,e_ratio(numz*numx*numz*numx,1), zt_value(numz*numx*numz*numx,1),zt_value2(numz*numx*numz*numx,1),x_value(numz*numx*numz*numx,1),omega_value(numz*numx*numz*numx,1),avg_merger_prob
    real*8, intent(out)::modelmom2(nmom,1)
    integer::errorflag

    real*8 :: decomp_diffv
    integer :: decomp_counter

    real*8::mean_z,dummy_z,std_z,mean_zt,dummy_zt,std_zt
    real*8::normalized_z(numz),normalized_zt(numz),normalized_zt2(numz)




    real*8,allocatable::initial_z_shock(:,:),psi_shock(:,:),eta_shock(:,:),innov_shock(:,:),innov_other_shock(:,:),z_other_shock(:,:),tiebreaking_shock(:,:),omega_shock(:,:), &
        merger_interp_shock(:,:), sim_equity_ratio(:,:)
    integer,allocatable::sim_owner(:,:),sim_z(:,:),sim_z_preinnov(:,:),sim_z_postinnov(:,:),sim_was_acquired_dummy(:,:),sim_rec_offer_dummy(:,:),sim_has_acquired_dummy(:,:)

    integer::   is_acquirer,imperfect_info

    real*8:: mu_lowerbar,pat_dep,mean_log_patstock,std_log_patstock
    real*8,allocatable:: innov(:,:),innov_back_dp(:,:),innov_stock(:,:),innov_stock2(:,:),innov_stock3(:,:),norm_log_patstock(:,:),norm_log_patstock_sq(:,:),bbeta_long(:,:),bbeta_long2(:,:),bbeta_long3(:,:),acq_dummy(:,:),innov_back_dp_vector(:,:)
    integer,allocatable::output_for_stata(:,:)
    real*8,allocatable::output_for_stata_real(:,:)
    integer:: regdim,regdim2,obs
    real*8 :: start, finish
    !double precision function omp_get_wtime()







    equity_ratio=0d0 
    equity_ratio_perfect=0d0 
    indexs_r = 0d0 
    indexc_r = 0d0 
    gc_r = 0d0 
    equity_ratio_r = 0d0 
    gacp_r = 0d0 
    merger_prob =0d0 
    conditional_accept = 0d0 
    equity_ratio_2 = 0d0 
    equity_ratio_3 = 0d0 
    offer_dummy = 0d0
    offer_dummy_2 = 0d0 
    corner_cash =0d0 
    corner_equity =0d0 
    interior =0d0 
    nonoffer =0d0 




    do i=1,numz  
        do j=1,numx 
            do ii=1,numz 
                do jj=1,numx

                    ! IMPERFECT INFORMATION CASE

                    inn = i_policy(ii,jj)
                    ! High type
                    iip = min(ii+1,numz)
                    x_vp = ((zt(iip)/zt(i))**phi * (zt(iip)+zt(i))**phi2)
                    g1s=LAMBDA1*( 1d0-gc_i(i,j,ii,jj,inn)/(gc_i(i,j,ii,jj,inn)+gs_i(i,j,ii,jj,inn)*(   max(0d0,pi*f(i,iip)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iip,j)*fv)      )) )**gammae2*f(i,iip) 
                    g1s=max(g1s,0d0)
                    g2s=LAMBDA1*( 1d0-gc_i(i,j,ii,jj,inn)/(gc_i(i,j,ii,jj,inn)+gs_i(i,j,ii,jj,inn)*(     max(0d0,pi*f(max(i-1,1),iip)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iip,j)*fv)       )) )**gammae2*f(max(i-1,1),iip) 
                    g2s=max(g2s,0d0)

                    temp1h = (pi*f(i,iip)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,iip,j))
                    temp1l = (pi*f(max(i-1,1),iip)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iip,j))
                    temp1h = max(temp1h,1e-10)
                    temp1l = max(temp1l,1e-10)
                    temp1s_h = kesi*temp1h + (1d0-kesi)*temp1l

                    ! Medium type
                    iip = ii
                    x_vp = ((zt(iip)/zt(i))**phi * (zt(iip)+zt(i))**phi2)
                    g1s=LAMBDA1*( 1d0-gc_i(i,j,ii,jj,inn)/(gc_i(i,j,ii,jj,inn)+gs_i(i,j,ii,jj,inn)*(   max(0d0,pi*f(i,iip)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iip,j)*fv)      )) )**gammae2*f(i,iip)
                    g1s=max(g1s,0d0)
                    g2s=LAMBDA1*( 1d0-gc_i(i,j,ii,jj,inn)/(gc_i(i,j,ii,jj,inn)+gs_i(i,j,ii,jj,inn)*(     max(0d0,pi*f(max(i-1,1),iip)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iip,j)*fv)       )) )**gammae2*f(max(i-1,1),iip) 
                    g2s=max(g2s,0d0)

                    temp1h = (pi*f(i,iip)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,iip,j))
                    temp1l = (pi*f(max(i-1,1),iip)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iip,j))
                    temp1h = max(temp1h,1e-10)
                    temp1l = max(temp1l,1e-10)
                    temp1s_m = kesi*temp1h + (1d0-kesi)*temp1l

                    ! Low type

                    iip = max(ii-1,1)
                    x_vp = ((zt(iip)/zt(i))**phi * (zt(iip)+zt(i))**phi2)
                    g1s=LAMBDA1*( 1d0-gc_i(i,j,ii,jj,inn)/(gc_i(i,j,ii,jj,inn)+gs_i(i,j,ii,jj,inn)*(   max(0d0,pi*f(i,iip)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iip,j)*fv)      )) )**gammae2*f(i,iip)
                    g1s=max(g1s,0d0)
                    g2s=LAMBDA1*( 1d0-gc_i(i,j,ii,jj,inn)/(gc_i(i,j,ii,jj,inn)+gs_i(i,j,ii,jj,inn)*(     max(0d0,pi*f(max(i-1,1),iip)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iip,j)*fv)       )) )**gammae2*f(max(i-1,1),iip) 
                    g2s=max(g2s,0d0)

                    temp1h = (pi*f(i,iip)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,iip,j))
                    temp1l = (pi*f(max(i-1,1),iip)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iip,j))
                    temp1h = max(temp1h,1e-10)
                    temp1l = max(temp1l,1e-10)
                    temp1s_l = kesi*temp1h + (1d0-kesi)*temp1l

                    temp1s = p3(1,ii,inn)*temp1s_h + p3(2,ii,inn)*temp1s_m + p3(3,ii,inn)*temp1s_l

                    equity_ratio(i,j,ii,jj)=1d0-gc_i(i,j,ii,jj,inn)/(gc_i(i,j,ii,jj,inn)+gs_i(i,j,ii,jj,inn)*temp1s)

                    ! PERFECT INFORMATION CASE

                    iip = ii
                    x_vp = ((zt(iip)/zt(i))**phi * (zt(iip)+zt(i))**phi2)
                    g1s=LAMBDA1*( 1d0-gc_perfect(i,j,ii,jj)/(gc_perfect(i,j,ii,jj)+gs_perfect(i,j,ii,jj)*(   max(0d0,pi*f(i,iip)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iip,j)*fv)      )) )**gammae2*f(i,iip)
                    g1s=max(g1s,0d0)
                    g2s=LAMBDA1*( 1d0-gc_perfect(i,j,ii,jj)/(gc_perfect(i,j,ii,jj)+gs_perfect(i,j,ii,jj)*(     max(0d0,pi*f(max(i-1,1),iip)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iip,j)*fv)       )) )**gammae2*f(max(i-1,1),iip) 
                    g2s=max(g2s,0d0)

                    temp1h = (pi*f(i,iip)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,iip,j))
                    temp1l = (pi*f(max(i-1,1),iip)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iip,j))
                    temp1h = max(temp1h,1e-10)
                    temp1l = max(temp1l,1e-10)
                    temp1s = kesi*temp1h + (1d0-kesi)*temp1l

                    equity_ratio_perfect(i,j,ii,jj)=1d0-gc_perfect(i,j,ii,jj)/(gc_perfect(i,j,ii,jj)+gs_perfect(i,j,ii,jj)*temp1s)

                end do
            end do
        end do
    end do

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



    helper_var1=0d0
    helper_var2=0d0
    helper_var3=0d0
    helper_equity=0d0
    new_equity_ratio_2=0d0
    new_corner_cash=0d0
    new_corner_equity=0d0
    new_interior=0d0
    new_nonoffer=0d0
    realized_gain=0d0
    realized_gain_helper=0d0
    acq_made_offer_helper=0d0

    realized_gain_dollar=0d0
    realized_gain_det=0d0
    realized_gain_det_helper=0d0
    equity_ratio_acq = 0d0
    value_loss = 0d0
    value_loss2 = 0d0
    acq_ann_ret = 0d0
    tar_ann_ret = 0d0
    tar_frac_merger_gain = 0d0
    tar_frac_of_gain = 0d0
    cash_ratio = 0d0
    log_mv_ratio = 0d0
    acq_ann_ret_helper = 0d0
    fixedcost_ratio=0d0
    rel_size = 0d0

    mve_reg_begin = 0d0
    mve_reg_end = 0d0
    mve_reg_helper = 0d0
    z_alone_begin = 0d0
    z_alone_end = 0d0
    z_alone_helper = 0d0

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NEW TRANS CODE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    trans = 0d0
    dummy = 0d0

    do j = 1,numx ! x
        do i = 1,numz ! z
            ! exogenous exit
            do ii_e = 1,numz ! entrant z
                do jj_e = 1,numx ! entrant x
                    call giveindex(i,j)
                    a=zx_idx !today index
                    call giveindex(ii_e,jj_e)
                    b=zx_idx !tomorrow index

                    F_exog_corrected = F_exog(ii_e,jj_e)/sum(F_exog)
                    if (sum(F_exog)==0d0) then
                        F_exog_corrected = 0d0
                    end if

                    trans(a,b) = trans(a,b) + psi*F_exog_corrected
                    dummy(i,j) = dummy(i,j) +  psi*F_exog_corrected
                end do
            end do
            ! no exit
            ! no meeting
            ! high z'
            call giveindex(i,j)
            a=zx_idx !today index
            call giveindex(min(i+1,numz),j)
            b=zx_idx !tomorrow index
            trans(a,b) = trans(a,b) + (1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
            dummy(i,j) = dummy(i,j) + (1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
            mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*(1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
            mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(min(i+1,numz),j)*(1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
            mve_reg_helper(i,j) = mve_reg_helper(i,j) + (1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
            z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*(1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
            z_alone_end(i,j) = z_alone_end(i,j) +  za(min(i+1,numz))*(1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
            z_alone_helper(i,j) = z_alone_helper(i,j) + (1d0-eta)*p3(1,i,i_policy(i,j))*(1d0-psi)
            ! med z'
            call giveindex(i,j)
            a=zx_idx !today index
            call giveindex(i,j)
            b=zx_idx !tomorrow index
            trans(a,b) = trans(a,b) + (1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
            dummy(i,j) = dummy(i,j) + (1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
            mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*(1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
            mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(i,j)*(1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
            mve_reg_helper(i,j) = mve_reg_helper(i,j) + (1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
            z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*(1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
            z_alone_end(i,j) = z_alone_end(i,j) + za(i)*(1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
            z_alone_helper(i,j) = z_alone_helper(i,j) + (1d0-eta)*p3(2,i,i_policy(i,j))*(1d0-psi)
            ! low z'
            call giveindex(i,j)
            a=zx_idx !today index
            call giveindex(max(i-1,1),j)
            b=zx_idx !tomorrow index
            trans(a,b) = trans(a,b) + (1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)
            dummy(i,j) = dummy(i,j) + (1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)
            mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*(1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)
            mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(max(i-1,1),j)*(1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)
            mve_reg_helper(i,j) = mve_reg_helper(i,j) + (1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)
            z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*(1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)
            z_alone_end(i,j) = z_alone_end(i,j) + za(max(i-1,1))*(1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)
            z_alone_helper(i,j) = z_alone_helper(i,j) + (1d0-eta)*p3(3,i,i_policy(i,j))*(1d0-psi)

            !if (be_acq(i,j) == 1) then ! Acquirer
            do ii = 1,numz ! other pre-innovation productivity
                do jj = 1,numx ! other fixed cost

                    if (i>ii) then
                        prob_acq(i,j) = 1d0
                    elseif (i==ii) then
                        prob_acq(i,j) = 0.5
                    else
                        prob_acq(i,j) = 0d0
                    endif

                    ! meet, acquirer -- imperfect information case
                    do k = -1,1 ! acquirer innovation realization
                        do m = -1,1 ! target innovation realization

                            ip = max(min(i+k,numz),1)
                            iip = max(min(ii+m,numz),1)
                            inn_a = i_policy(i,j)
                            inn_t = i_policy(ii,jj)
                            tttr = Fs(ii,jj)/sum(Fs)*info_asym*prob_acq(i,j) *eta*p3(2-k,i,inn_a)*p3(2-m,ii,inn_t)*(1d0-psi)

                            x_vp = ((zt(iip)/zt(ip))**phi * (zt(iip)+zt(ip))**phi2)


                            if (sum(Fs)==0d0) then
                                tttr = 0d0
                            end if
                            if (offermade_i(ip,j,ii,jj,i_policy(ii,jj)) == 1) then


                                acq_ann_ret(i,j,ii,jj) = acq_ann_ret(i,j,ii,jj) + ((Vanew_i(ip,j,ii,jj,inn_t)-Vnew(i,j))/Vnew(i,j))*tttr
                                if (m == 1) then
                                    temp_tar_ann_ret = (Vtnew_h_i(ip,j,iip,jj,inn_t)-Vnew(ii,jj))/(Vnew(ii,jj))
                                elseif (m == 0) then
                                    temp_tar_ann_ret = (Vtnew_m_i(ip,j,iip,jj,inn_t)-Vnew(ii,jj))/(Vnew(ii,jj))
                                else
                                    temp_tar_ann_ret = (Vtnew_l_i(ip,j,iip,jj,inn_t)-Vnew(ii,jj))/(Vnew(ii,jj))
                                endif
                                tar_ann_ret(i,j,ii,jj) = tar_ann_ret(i,j,ii,jj) + temp_tar_ann_ret*tttr
                                tar_frac_merger_gain(i,j,ii,jj) = tar_frac_merger_gain(i,j,ii,jj) + ( (  (temp_tar_ann_ret*Vnew(ii,jj)) - (Vanew_i(ip,j,ii,jj,inn_t)-Vnew(i,j)) ) / ( Vnew(ii,jj) + Vnew(i,j)) )*tttr
                                cash_ratio(i,j,ii,jj) =  cash_ratio(i,j,ii,jj) +  (1d0-equity_ratio(ip,j,ii,jj))*tttr
                                log_mv_ratio(i,j,ii,jj) = log_mv_ratio(i,j,ii,jj) + (LOG(Vnew(i,j)/Vnew(ii,jj)))*tttr
                                rel_size(i,j,ii,jj) = rel_size(i,j,ii,jj) + (za(ii)/za(i))*tttr

                                acq_ann_ret_helper(i,j,ii,jj) = acq_ann_ret_helper(i,j,ii,jj) + tttr
                                acq_made_offer_helper(i,j) = acq_made_offer_helper(i,j) + tttr

                                if (gacp_eqm_i(ip,j,ii,jj,2-m,i_policy(ii,jj)) == 1) then! if target accepts
                                    ! merger risk -- no drop
                                    call zinterp(ip,iip,p_lb, p_ub,i_lb,i_ub)
                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(i_lb,j)
                                    b=zx_idx !tomorrow index

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*kesi*p_lb*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(i_lb,j)*kesi*p_lb*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + kesi*p_lb*tttr


                                    g1s=LAMBDA1*( 1d0-gc_i(ip,j,ii,jj,inn_t)/(gc_i(ip,j,ii,jj,inn_t)+gs_i(ip,j,ii,jj,inn_t)*(  max(0d0,pi*za(i_lb)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V(i_lb,j)*fv)      )) )**gammae2*za(i_lb)
                                    g1s=max(g1s,0d0)

                                    merged_company_value = (pi*za(i_lb) - x(j)*x_vp - g1s + (1d0-psi)/(1d0+r)*V(i_lb,j))
                                    premerge_acq_value = (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))
                                    premerge_tar_value = (pi*za(ii)+(1d0-psi)/(1d0+r)*V(ii,jj)) !based on observed type of the target
                                    realized_gain_current = merged_company_value / (premerge_acq_value+premerge_tar_value)
                                    realized_gain_current_dollar = merged_company_value - (premerge_acq_value+premerge_tar_value)
                                    tar_frac_of_gain(i,j) = tar_frac_of_gain(i,j) +  kesi*p_lb*tttr*((temp_tar_ann_ret*(Vnew(ii,jj)))+Vnew(ii,jj)-Vnew(ii,jj))/(merged_company_value-(V(i,j)+V(ii,jj)))


                                    fixedcost_ratio(i,j)=fixedcost_ratio(i,j)+  kesi*p_lb*tttr * (x(j)*x_vp/merged_company_value)
                                    realized_gain(i,j) = realized_gain(i,j) + kesi*p_lb*tttr * realized_gain_current
                                    realized_gain_dollar(i,j) = realized_gain_dollar(i,j) + kesi*p_lb*tttr * realized_gain_current_dollar

                                    realized_gain_helper(i,j) = realized_gain_helper(i,j) + kesi*p_lb*tttr


                                    realized_gain_det(ip,j,iip,jj,1,1) = realized_gain_det(ip,j,iip,jj,1,1) + kesi*p_lb*tttr * realized_gain_current
                                    realized_gain_det_helper(ip,j,iip,jj,1,1) = realized_gain_det_helper(ip,j,iip,jj,1,1) + kesi*p_lb*tttr

                                    equity_ratio_acq(ip,j,iip,jj,1,1) = equity_ratio_acq(ip,j,iip,jj,1,1) + kesi*p_lb*tttr * equity_ratio(ip,j,ii,jj)
                                    value_loss(ip,j,iip,jj,1,1) = value_loss(ip,j,iip,jj,1,1) + kesi*p_lb*tttr * (g1s/merged_company_value)
                                    value_loss2(ip,j,iip,jj,1,1) = value_loss2(ip,j,iip,jj,1,1) + kesi*p_lb*tttr * (g1s)

                                    trans(a,b)= trans(a,b)+ kesi*p_lb*tttr

                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(i_ub,j)
                                    b=zx_idx !tomorrow index

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*kesi*p_ub*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(i_ub,j)*kesi*p_ub*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + kesi*p_ub*tttr

                                    g2s=LAMBDA1*( 1d0-gc_i(ip,j,ii,jj,inn_t)/(gc_i(ip,j,ii,jj,inn_t)+gs_i(ip,j,ii,jj,inn_t)*(  max(0d0,pi*za(i_ub)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V(i_ub,j)*fv)     )) )**gammae2*za(i_ub)
                                    g2s=max(g2s,0d0)

                                    merged_company_value = (pi*za(i_ub) - x(j)*x_vp - g2s + (1d0-psi)/(1d0+r)*V(i_ub,j))
                                    premerge_acq_value = (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))
                                    premerge_tar_value = (pi*za(ii)+(1d0-psi)/(1d0+r)*V(ii,jj)) !based on observed type of the target
                                    realized_gain_current = merged_company_value / (premerge_acq_value+premerge_tar_value)
                                    realized_gain_current_dollar = merged_company_value - (premerge_acq_value+premerge_tar_value)
                                    tar_frac_of_gain(i,j) = tar_frac_of_gain(i,j) +  kesi*p_ub*tttr*((temp_tar_ann_ret*(Vnew(ii,jj)))+Vnew(ii,jj)-Vnew(ii,jj))/(merged_company_value-(V(i,j)+V(ii,jj)))


                                    fixedcost_ratio(i,j)=fixedcost_ratio(i,j)+   kesi*p_ub*tttr * (x(j)*x_vp/merged_company_value)
                                    realized_gain(i,j) = realized_gain(i,j) + kesi*p_ub*tttr * realized_gain_current
                                    realized_gain_dollar(i,j) = realized_gain_dollar(i,j) + kesi*p_ub*tttr * realized_gain_current_dollar

                                    realized_gain_helper(i,j) = realized_gain_helper(i,j) + kesi*p_ub*tttr

                                    realized_gain_det(ip,j,iip,jj,2,1) = realized_gain_det(ip,j,iip,jj,2,1) + kesi*p_ub*tttr * realized_gain_current
                                    realized_gain_det_helper(ip,j,iip,jj,2,1) = realized_gain_det_helper(ip,j,iip,jj,2,1) + kesi*p_ub*tttr

                                    equity_ratio_acq(ip,j,iip,jj,2,1) = equity_ratio_acq(ip,j,iip,jj,2,1) + kesi*p_ub*tttr * equity_ratio(ip,j,ii,jj)
                                    value_loss(ip,j,iip,jj,2,1) = value_loss(ip,j,iip,jj,2,1) + kesi*p_ub*tttr * (g2s/merged_company_value)
                                    value_loss2(ip,j,iip,jj,2,1) = value_loss2(ip,j,iip,jj,2,1) + kesi*p_ub*tttr * (g2s)


                                    trans(a,b)= trans(a,b)+ kesi*p_ub*tttr
                                    dummy(i,j) = dummy(i,j) + kesi*p_lb*tttr
                                    dummy(i,j) = dummy(i,j) + kesi*p_ub*tttr

                                    ! merger risk -- drop
                                    !   [ p_lb, p_ub,i_lb,i_ub]=zinterp(max(ip-1,1),iip)

                                    call zinterp(max(ip-1,1),iip,p_lb, p_ub,i_lb,i_ub)
                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(i_lb,j)
                                    b=zx_idx !tomorrow index

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*(1d0-kesi)*p_lb*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(i_lb,j)*(1d0-kesi)*p_lb*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + (1d0-kesi)*p_lb*tttr

                                    g1s=LAMBDA1*( 1d0-gc_i(ip,j,ii,jj,inn_t)/(gc_i(ip,j,ii,jj,inn_t)+gs_i(ip,j,ii,jj,inn_t)*(  max(0d0,pi*za(i_lb)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V(i_lb,j)*fv)      )) )**gammae2*za(i_lb)
                                    g1s=max(g1s,0d0)

                                    merged_company_value = (pi*za(i_lb) - x(j)*x_vp - g1s + (1d0-psi)/(1d0+r)*V(i_lb,j))
                                    premerge_acq_value = (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))
                                    premerge_tar_value = (pi*za(ii)+(1d0-psi)/(1d0+r)*V(ii,jj)) !based on observed type of the target
                                    realized_gain_current = merged_company_value / (premerge_acq_value+premerge_tar_value)
                                    realized_gain_current_dollar = merged_company_value - (premerge_acq_value+premerge_tar_value)

                                    tar_frac_of_gain(i,j) = tar_frac_of_gain(i,j) +  (1d0-kesi)*p_lb*tttr*((temp_tar_ann_ret*(Vnew(ii,jj)))+Vnew(ii,jj)-Vnew(ii,jj))/(merged_company_value-(V(i,j)+V(ii,jj)))


                                    fixedcost_ratio(i,j)=fixedcost_ratio(i,j)+   (1d0-kesi)*p_lb*tttr *  (x(j)*x_vp/merged_company_value)
                                    realized_gain(i,j) = realized_gain(i,j) + (1d0-kesi)*p_lb*tttr * realized_gain_current
                                    realized_gain_dollar(i,j) = realized_gain_dollar(i,j) + (1d0-kesi)*p_lb*tttr * realized_gain_current_dollar

                                    realized_gain_helper(i,j) = realized_gain_helper(i,j) + (1d0-kesi)*p_lb*tttr

                                    realized_gain_det(ip,j,iip,jj,1,2) = realized_gain_det(ip,j,iip,jj,1,2) +  (1d0-kesi)*p_lb*tttr * realized_gain_current
                                    realized_gain_det_helper(ip,j,iip,jj,1,2) = realized_gain_det_helper(ip,j,iip,jj,1,2) + (1d0-kesi)*p_lb*tttr

                                    equity_ratio_acq(ip,j,iip,jj,1,2) = equity_ratio_acq(ip,j,iip,jj,1,2) + (1d0-kesi)*p_lb*tttr * equity_ratio(ip,j,ii,jj)
                                    value_loss(ip,j,iip,jj,1,2) = value_loss(ip,j,iip,jj,1,2) + (1d0-kesi)*p_lb*tttr * (g1s/merged_company_value)
                                    value_loss2(ip,j,iip,jj,1,2) = value_loss2(ip,j,iip,jj,1,2) + (1d0-kesi)*p_lb*tttr * (g1s)


                                    trans(a,b)= trans(a,b)+ (1d0-kesi)*p_lb*tttr
                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(i_ub,j)
                                    b=zx_idx !tomorrow index

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*(1d0-kesi)*p_ub*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(i_ub,j)*(1d0-kesi)*p_ub*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + (1d0-kesi)*p_ub*tttr


                                    g2s=LAMBDA1*( 1d0-gc_i(ip,j,ii,jj,inn_t)/(gc_i(ip,j,ii,jj,inn_t)+gs_i(ip,j,ii,jj,inn_t)*(  max(0d0,pi*za(i_ub)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V(i_ub,j)*fv)     )) )**gammae2*za(i_ub)
                                    g2s=max(g2s,0d0)

                                    merged_company_value = (pi*za(i_ub) - x(j)*x_vp - g2s + (1d0-psi)/(1d0+r)*V(i_ub,j))
                                    premerge_acq_value = (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))
                                    premerge_tar_value = (pi*za(ii)+(1d0-psi)/(1d0+r)*V(ii,jj)) !based on observed type of the target
                                    realized_gain_current = merged_company_value / (premerge_acq_value+premerge_tar_value)
                                    realized_gain_current_dollar = merged_company_value - (premerge_acq_value+premerge_tar_value)

                                    tar_frac_of_gain(i,j) = tar_frac_of_gain(i,j) +  (1d0-kesi)*p_ub*tttr*((temp_tar_ann_ret*(Vnew(ii,jj)))+Vnew(ii,jj)-Vnew(ii,jj))/(merged_company_value-(V(i,j)+V(ii,jj)))


                                    fixedcost_ratio(i,j)=fixedcost_ratio(i,j)+   (1d0-kesi)*p_ub*tttr * (x(j)*x_vp/merged_company_value)
                                    realized_gain(i,j) = realized_gain(i,j) + (1d0-kesi)*p_ub*tttr * realized_gain_current
                                    realized_gain_dollar(i,j) = realized_gain_dollar(i,j) + (1d0-kesi)*p_ub*tttr * realized_gain_current_dollar

                                    realized_gain_helper(i,j) = realized_gain_helper(i,j) + (1d0-kesi)*p_ub*tttr

                                    realized_gain_det(ip,j,iip,jj,2,2) = realized_gain_det(ip,j,iip,jj,2,2) +  (1d0-kesi)*p_ub*tttr * realized_gain_current
                                    realized_gain_det_helper(ip,j,iip,jj,2,2) = realized_gain_det_helper(ip,j,iip,jj,2,2) + (1d0-kesi)*p_ub*tttr

                                    equity_ratio_acq(ip,j,iip,jj,2,2) = equity_ratio_acq(ip,j,iip,jj,2,2) + (1d0-kesi)*p_ub*tttr * equity_ratio(ip,j,ii,jj)
                                    value_loss(ip,j,iip,jj,2,2) = value_loss(ip,j,iip,jj,2,2) + (1d0-kesi)*p_ub*tttr * (g2s/merged_company_value)
                                    value_loss2(ip,j,iip,jj,2,2) = value_loss2(ip,j,iip,jj,2,2) + (1d0-kesi)*p_ub*tttr * (g2s)


                                    trans(a,b)= trans(a,b)+ (1d0-kesi)*p_ub*tttr
                                    dummy(i,j) = dummy(i,j) + (1d0-kesi)*p_lb*tttr
                                    dummy(i,j) = dummy(i,j) + (1d0-kesi)*p_ub*tttr

                                else ! if target rejects


                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call  giveindex(ip,j)
                                    b=zx_idx !tomorrow index


                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(ip,j)*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + tttr
                                    z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*tttr
                                    z_alone_end(i,j) = z_alone_end(i,j) + za(ip)*tttr
                                    z_alone_helper(i,j) = z_alone_helper(i,j) + tttr


                                    trans(a,b) = trans(a,b) + tttr
                                    dummy(i,j) = dummy(i,j) +  tttr
                                end if
                            else ! no offer
                                call giveindex(i,j)
                                a=zx_idx !today index
                                call  giveindex(ip,j)
                                b=zx_idx !tomorrow index

                                mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*tttr
                                mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(ip,j)*tttr
                                mve_reg_helper(i,j) = mve_reg_helper(i,j) + tttr
                                z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*tttr
                                z_alone_end(i,j) = z_alone_end(i,j) + za(ip)*tttr
                                z_alone_helper(i,j) = z_alone_helper(i,j) + tttr


                                trans(a,b) = trans(a,b) + tttr
                                dummy(i,j) = dummy(i,j) + tttr
                            end if
                        end do
                    end do
                end do
            end do

            ! meet -- perfect information case
            do ii = 1,numz ! other pre-innovation productivity
                do jj = 1,numx ! other fixed cost

                    if (i>ii) then
                        prob_acq(i,j) = 1d0
                    elseif (i==ii) then
                        prob_acq(i,j) = 0.5
                    else
                        prob_acq(i,j) = 0d0
                    endif

                    do k = -1,1 ! acquirer innovation realization
                        do m = -1,1 ! target innovation realization
                            ip = max(min(i+k,numz),1)
                            iip = max(min(ii+m,numz),1)

                            inn_a = i_policy(i,j)
                            inn_t = i_policy(ii,jj)
                            tttr = Fs(ii,jj)/sum(Fs)*(1d0-info_asym)*prob_acq(i,j) * eta*p3(2-k,i,inn_a)*p3(2-m,ii,inn_t)*(1d0-psi)

                            x_vp = ((zt(iip)/zt(ip))**phi * (zt(iip)+zt(ip))**phi2)

                            if (sum(Fs)==0d0) then
                                tttr = 0d0
                            end if



                            if (offermade_perfect(ip,j,iip,jj) == 1) then

                                acq_ann_ret(i,j,iip,jj) = acq_ann_ret(i,j,iip,jj) + ((Vanew_perfect(ip,j,iip,jj)-Vnew(i,j))/Vnew(i,j))*tttr
                                temp_tar_ann_ret = (Vtnew_perfect(ip,j,iip,jj)-Vnew(ii,jj))/(Vnew(ii,jj))
                                tar_ann_ret(i,j,ii,jj) = tar_ann_ret(i,j,ii,jj) + temp_tar_ann_ret*tttr
                                tar_frac_merger_gain(i,j,ii,jj) = tar_frac_merger_gain(i,j,ii,jj) + ( (  (temp_tar_ann_ret*Vnew(ii,jj)) - (Vanew_perfect(ip,j,iip,jj)-Vnew(i,j)) ) / ( Vnew(ii,jj) + Vnew(i,j)) )*tttr
                                cash_ratio(i,j,iip,jj) =  cash_ratio(i,j,iip,jj) +  (1d0-equity_ratio_perfect(ip,j,iip,jj))*tttr
                                log_mv_ratio(i,j,iip,jj) = log_mv_ratio(i,j,iip,jj) + (LOG(Vnew(i,j)/Vnew(ii,jj)))*tttr
                                rel_size(i,j,ii,jj) = rel_size(i,j,ii,jj) + (za(ii)/za(i))*tttr


                                acq_ann_ret_helper(i,j,iip,jj) = acq_ann_ret_helper(i,j,iip,jj) + tttr
                                acq_made_offer_helper(i,j) = acq_made_offer_helper(i,j) + tttr

                                if (gacp_eqm_perfect(ip,j,iip,jj) == 1) then! if target accepts
                                    ! merger risk -- no drop
                                    call zinterp(ip,iip,p_lb, p_ub,i_lb,i_ub)
                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(i_lb,j)
                                    b=zx_idx !tomorrow index

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*kesi*p_lb*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(i_lb,j)*kesi*p_lb*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + kesi*p_lb*tttr

                                    g1s=LAMBDA1*( 1d0-gc_perfect(ip,j,iip,jj)/(gc_perfect(ip,j,iip,jj)+gs_perfect(ip,j,iip,jj)*(  max(0d0,pi*za(i_lb)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V(i_lb,j)*fv)      )) )**gammae2*za(i_lb)
                                    g1s=max(g1s,0d0)

                                    merged_company_value = (pi*za(i_lb) - x(j)*x_vp - g1s + (1d0-psi)/(1d0+r)*V(i_lb,j))
                                    premerge_acq_value = (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))
                                    premerge_tar_value = (pi*za(iip)+(1d0-psi)/(1d0+r)*V(iip,jj)) 
                                    realized_gain_current = merged_company_value / (premerge_acq_value+premerge_tar_value)
                                    realized_gain_current_dollar = merged_company_value - (premerge_acq_value+premerge_tar_value)

                                    tar_frac_of_gain(i,j) = tar_frac_of_gain(i,j) +  kesi*p_lb*tttr*((temp_tar_ann_ret*(Vnew(ii,jj)))+Vnew(ii,jj)-Vnew(ii,jj))/(merged_company_value-(V(i,j)+V(ii,jj)))


                                    fixedcost_ratio(i,j)=fixedcost_ratio(i,j)+   kesi*p_lb*tttr * (x(j)*x_vp/merged_company_value)
                                    realized_gain(i,j) = realized_gain(i,j) + kesi*p_lb*tttr * realized_gain_current
                                    realized_gain_dollar(i,j) = realized_gain_dollar(i,j) + kesi*p_lb*tttr * realized_gain_current_dollar

                                    realized_gain_helper(i,j) = realized_gain_helper(i,j) + kesi*p_lb*tttr

                                    ! argument 5 = 1 means p_lb, else p_ub
                                    ! argument 6 = 1 means kesi, else 1-kesi
                                    realized_gain_det(ip,j,iip,jj,1,1) = realized_gain_det(ip,j,iip,jj,1,1) + kesi*p_lb*tttr * realized_gain_current
                                    realized_gain_det_helper(ip,j,iip,jj,1,1) = realized_gain_det_helper(ip,j,iip,jj,1,1) + kesi*p_lb*tttr

                                    equity_ratio_acq(ip,j,iip,jj,1,1) = equity_ratio_acq(ip,j,iip,jj,1,1) + kesi*p_lb*tttr * equity_ratio_perfect(ip,j,iip,jj)
                                    value_loss(ip,j,iip,jj,1,1) = value_loss(ip,j,iip,jj,1,1) + kesi*p_lb*tttr * (g1s/merged_company_value)
                                    value_loss2(ip,j,iip,jj,1,1) = value_loss2(ip,j,iip,jj,1,1) + kesi*p_lb*tttr * (g1s)


                                    trans(a,b)= trans(a,b)+ kesi*p_lb*tttr


                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(i_ub,j)
                                    b=zx_idx !tomorrow index

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*kesi*p_ub*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(i_ub,j)*kesi*p_ub*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + kesi*p_ub*tttr

                                    g2s=LAMBDA1*( 1d0-gc_perfect(ip,j,iip,jj)/(gc_perfect(ip,j,iip,jj)+gs_perfect(ip,j,iip,jj)*(  max(0d0,pi*za(i_ub)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V(i_ub,j)*fv)     )) )**gammae2*za(i_ub)
                                    g2s=max(g2s,0d0)

                                    merged_company_value = (pi*za(i_ub) - x(j)*x_vp - g2s + (1d0-psi)/(1d0+r)*V(i_ub,j))
                                    premerge_acq_value = (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))
                                    premerge_tar_value = (pi*za(iip)+(1d0-psi)/(1d0+r)*V(iip,jj)) 
                                    realized_gain_current = merged_company_value / (premerge_acq_value+premerge_tar_value)
                                    realized_gain_current_dollar = merged_company_value - (premerge_acq_value+premerge_tar_value)

                                    tar_frac_of_gain(i,j) = tar_frac_of_gain(i,j) +  kesi*p_ub*tttr*((temp_tar_ann_ret*(Vnew(ii,jj)))+Vnew(ii,jj)-Vnew(ii,jj))/(merged_company_value-(V(i,j)+V(ii,jj)))


                                    fixedcost_ratio(i,j)=fixedcost_ratio(i,j)+ kesi*p_ub*tttr *  (x(j)*x_vp/merged_company_value)
                                    realized_gain(i,j) = realized_gain(i,j) + kesi*p_ub*tttr * realized_gain_current
                                    realized_gain_dollar(i,j) = realized_gain_dollar(i,j) + kesi*p_ub*tttr * realized_gain_current_dollar

                                    realized_gain_helper(i,j) = realized_gain_helper(i,j) + kesi*p_ub*tttr

                                    realized_gain_det(ip,j,iip,jj,2,1) = realized_gain_det(ip,j,iip,jj,2,1) + kesi*p_ub*tttr * realized_gain_current
                                    realized_gain_det_helper(ip,j,iip,jj,2,1) = realized_gain_det_helper(ip,j,iip,jj,2,1) + kesi*p_ub*tttr

                                    equity_ratio_acq(ip,j,iip,jj,2,1) = equity_ratio_acq(ip,j,iip,jj,2,1) + kesi*p_ub*tttr * equity_ratio_perfect(ip,j,iip,jj)
                                    value_loss(ip,j,iip,jj,2,1) = value_loss(ip,j,iip,jj,2,1) + kesi*p_ub*tttr * (g2s/merged_company_value)
                                    value_loss2(ip,j,iip,jj,2,1) = value_loss2(ip,j,iip,jj,2,1) + kesi*p_ub*tttr * (g2s)


                                    trans(a,b)= trans(a,b)+ kesi*p_ub*tttr
                                    dummy(i,j) = dummy(i,j) + kesi*p_lb*tttr
                                    dummy(i,j) = dummy(i,j) + kesi*p_ub*tttr

                                    ! merger risk -- drop
                                    !   [ p_lb, p_ub,i_lb,i_ub]=zinterp(max(ip-1,1),iip)

                                    call zinterp(max(ip-1,1),iip,p_lb, p_ub,i_lb,i_ub)
                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(i_lb,j)
                                    b=zx_idx !tomorrow index

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*(1d0-kesi)*p_lb*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(i_lb,j)*(1d0-kesi)*p_lb*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + (1d0-kesi)*p_lb*tttr

                                    g1s=LAMBDA1*( 1d0-gc_perfect(ip,j,iip,jj)/(gc_perfect(ip,j,iip,jj)+gs_perfect(ip,j,iip,jj)*(  max(0d0,pi*za(i_lb)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V(i_lb,j)*fv)      )) )**gammae2*za(i_lb)
                                    g1s=max(g1s,0d0)

                                    merged_company_value = (pi*za(i_lb) - x(j)*x_vp - g1s + (1d0-psi)/(1d0+r)*V(i_lb,j))
                                    premerge_acq_value = (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))
                                    premerge_tar_value = (pi*za(iip)+(1d0-psi)/(1d0+r)*V(iip,jj)) 
                                    realized_gain_current = merged_company_value / (premerge_acq_value+premerge_tar_value)
                                    realized_gain_current_dollar = merged_company_value - (premerge_acq_value+premerge_tar_value)

                                    tar_frac_of_gain(i,j) = tar_frac_of_gain(i,j) +  (1d0-kesi)*p_lb*tttr*((temp_tar_ann_ret*(Vnew(ii,jj)))+Vnew(ii,jj)-Vnew(ii,jj))/(merged_company_value-(V(i,j)+V(ii,jj)))

                                    fixedcost_ratio(i,j)=fixedcost_ratio(i,j)+   (1d0-kesi)*p_lb*tttr * (x(j)*x_vp/merged_company_value)
                                    realized_gain(i,j) = realized_gain(i,j) + (1d0-kesi)*p_lb*tttr * realized_gain_current
                                    realized_gain_dollar(i,j) = realized_gain_dollar(i,j) + (1d0-kesi)*p_lb*tttr * realized_gain_current_dollar

                                    realized_gain_helper(i,j) = realized_gain_helper(i,j) + (1d0-kesi)*p_lb*tttr

                                    realized_gain_det(ip,j,iip,jj,1,2) = realized_gain_det(ip,j,iip,jj,1,2) +  (1d0-kesi)*p_lb*tttr * realized_gain_current
                                    realized_gain_det_helper(ip,j,iip,jj,1,2) = realized_gain_det_helper(ip,j,iip,jj,1,2) + (1d0-kesi)*p_lb*tttr

                                    equity_ratio_acq(ip,j,iip,jj,1,2) = equity_ratio_acq(ip,j,iip,jj,1,2) +(1d0-kesi)*p_lb*tttr * equity_ratio_perfect(ip,j,iip,jj)
                                    value_loss(ip,j,iip,jj,1,2) = value_loss(ip,j,iip,jj,1,2) + (1d0-kesi)*p_lb*tttr * (g1s/merged_company_value)
                                    value_loss2(ip,j,iip,jj,1,2) = value_loss2(ip,j,iip,jj,1,2) + (1d0-kesi)*p_lb*tttr * (g1s)


                                    trans(a,b)= trans(a,b)+ (1d0-kesi)*p_lb*tttr
                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(i_ub,j)
                                    b=zx_idx !tomorrow index

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*(1d0-kesi)*p_ub*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(i_ub,j)*(1d0-kesi)*p_ub*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + (1d0-kesi)*p_ub*tttr

                                    g2s=LAMBDA1*( 1d0-gc_perfect(ip,j,iip,jj)/(gc_perfect(ip,j,iip,jj)+gs_perfect(ip,j,iip,jj)*(  max(0d0,pi*za(i_ub)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V(i_ub,j)*fv)     )) )**gammae2*za(i_ub)
                                    g2s=max(g2s,0d0)

                                    merged_company_value = (pi*za(i_ub) - x(j)*x_vp - g2s + (1d0-psi)/(1d0+r)*V(i_ub,j))
                                    premerge_acq_value = (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))
                                    premerge_tar_value = (pi*za(iip)+(1d0-psi)/(1d0+r)*V(iip,jj)) 
                                    realized_gain_current = merged_company_value / (premerge_acq_value+premerge_tar_value)
                                    realized_gain_current_dollar = merged_company_value - (premerge_acq_value+premerge_tar_value)

                                    tar_frac_of_gain(i,j) = tar_frac_of_gain(i,j) +  (1d0-kesi)*p_ub*tttr*((temp_tar_ann_ret*(Vnew(ii,jj)))+Vnew(ii,jj)-Vnew(ii,jj))/(merged_company_value-(V(i,j)+V(ii,jj)))

                                    fixedcost_ratio(i,j)=fixedcost_ratio(i,j)+  (1d0-kesi)*p_ub*tttr * (x(j)*x_vp/merged_company_value)
                                    realized_gain(i,j) = realized_gain(i,j) + (1d0-kesi)*p_ub*tttr * realized_gain_current
                                    realized_gain_dollar(i,j) = realized_gain_dollar(i,j) + (1d0-kesi)*p_ub*tttr * realized_gain_current_dollar

                                    realized_gain_helper(i,j) = realized_gain_helper(i,j) + (1d0-kesi)*p_ub*tttr

                                    realized_gain_det(ip,j,iip,jj,2,2) = realized_gain_det(ip,j,iip,jj,2,2) +  (1d0-kesi)*p_ub*tttr * realized_gain_current
                                    realized_gain_det_helper(ip,j,iip,jj,2,2) = realized_gain_det_helper(ip,j,iip,jj,2,2) + (1d0-kesi)*p_ub*tttr

                                    equity_ratio_acq(ip,j,iip,jj,2,2) = equity_ratio_acq(ip,j,iip,jj,2,2) + (1d0-kesi)*p_ub*tttr * equity_ratio_perfect(ip,j,iip,jj)
                                    value_loss(ip,j,iip,jj,2,2) = value_loss(ip,j,iip,jj,2,2) + (1d0-kesi)*p_ub*tttr * (g2s/merged_company_value)
                                    value_loss2(ip,j,iip,jj,2,2) = value_loss2(ip,j,iip,jj,2,2) + (1d0-kesi)*p_ub*tttr * (g2s)

                                    trans(a,b)= trans(a,b)+ (1d0-kesi)*p_ub*tttr
                                    dummy(i,j) = dummy(i,j) +  (1d0-kesi)*p_lb*tttr
                                    dummy(i,j) = dummy(i,j) + (1d0-kesi)*p_ub*tttr


                                else ! if target rejects


                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call  giveindex(ip,j)
                                    b=zx_idx !tomorrow index

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*tttr
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(ip,j)*tttr
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + tttr
                                    z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*tttr
                                    z_alone_end(i,j) = z_alone_end(i,j) + za(ip)*tttr
                                    z_alone_helper(i,j) = z_alone_helper(i,j) + tttr


                                    trans(a,b) = trans(a,b) + tttr
                                    dummy(i,j) = dummy(i,j) +  tttr
                                end if
                            else ! no offer
                                call giveindex(i,j)
                                a=zx_idx !today index
                                call  giveindex(ip,j)
                                b=zx_idx !tomorrow index

                                mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*tttr
                                mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(ip,j)*tttr
                                mve_reg_helper(i,j) = mve_reg_helper(i,j) + tttr
                                z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*tttr
                                z_alone_end(i,j) = z_alone_end(i,j) + za(ip)*tttr
                                z_alone_helper(i,j) = z_alone_helper(i,j) + tttr

                                trans(a,b) = trans(a,b) + tttr
                                dummy(i,j) = dummy(i,j) + tttr
                            end if
                        end do
                    end do
                end do
            end do

            !else ! Target
            ! meet -- imperfect information case
            do ii = 1,numz ! other pre-innovaton productivity
                do jj = 1,numx ! other fixed cost

                    if (i>ii) then
                        prob_acq(i,j) = 1d0
                    elseif (i==ii) then
                        prob_acq(i,j) = 0.5
                    else
                        prob_acq(i,j) = 0d0
                    endif

                    do k = -1,1 ! my innovation realization
                        do m = -1,1 ! other innovation realization

                            ip = max(min(i+k,numz),1) ! target after innov
                            iip = max(min(ii+m,numz),1) ! acquirer after innov

                            inn_a = i_policy(ii,jj)
                            inn_t = i_policy(i,j)
                            aaa = Fs(ii,jj)/sum(Fs)*info_asym*(1d0-prob_acq(i,j)) *eta*p3(2-k,i,inn_t)*p3(2-m,ii,inn_a)*(1d0-psi)
                            if (sum(Fs) == 0d0) then
                                aaa = 0d0
                            end if

                            ! helper_var2 records the probability of a meeting conditional on target z and x
                            helper_var2(i,j) = helper_var2(i,j) + aaa
                            if (offermade_i(iip,jj,i,j,i_policy(i,j)) == 1) then

                                ! helper_var3 records the probability of an offer being made conditional on target z and x and a meeting
                                helper_var3(i,j) = helper_var3(i,j) +aaa
                                helper_equity(i,j) = helper_equity(i,j) + aaa*equity_ratio(iip,jj,i,j)

                                if (indexs_i(iip,jj,i,j,inn_t)==1 .and. gc_i(iip,jj,i,j,inn_t)>1e-3) then
                                    new_corner_equity(i,j) = new_corner_equity(i,j) + aaa !;
                                end if

                                if (indexs_i(iip,jj,i,j,inn_t)/=1 .and. gc_i(iip,jj,i,j,inn_t)<=1e-3) then
                                    new_corner_cash(i,j) = new_corner_cash(i,j) + aaa !;
                                end if

                                if (indexs_i(iip,jj,i,j,inn_t)/=1 .and. gc_i(iip,jj,i,j,inn_t)>1e-3) then
                                    new_interior(i,j) = new_interior(i,j) + aaa !;
                                end if

                                if (indexs_i(iip,jj,i,j,inn_t)==1 .and. gc_i(iip,jj,i,j,inn_t)<=1e-3) then
                                    new_nonoffer(i,j) = new_nonoffer(i,j) + aaa!;
                                end if

                                if (gacp_eqm_i(iip,jj,i,j,2-k,i_policy(i,j)) == 1) then ! if target accepts


                                    ! helper_var1 records the probability of acceptance conditional on target z and x, a meeting, and an offer made

                                    helper_var1(i,j) = helper_var1(i,j) + aaa
                                    helper_equity1(i,j) = helper_equity1(i,j) +   aaa*equity_ratio(iip,jj,i,j)


                                    do ii_e = 1,numz ! entrant z
                                        do jj_e = 1,numx ! entrant x
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(ii_e,jj_e)
                                            b=zx_idx

                                            F_exog_corrected = F_exog(ii_e,jj_e)/sum(F_exog)
                                            if (sum(F_exog)==0d0) then
                                                F_exog_corrected = 0d0
                                            end if


                                            trans(a,b) = trans(a,b) +  aaa*(F_exog_corrected)
                                            dummy(i,j) = dummy(i,j) +  aaa*(F_exog_corrected)
                                        end do
                                    end do
                                else ! if target rejects
                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(ip,j)
                                    b=zx_idx

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*aaa
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(ip,j)*aaa
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + aaa
                                    z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*aaa
                                    z_alone_end(i,j) = z_alone_end(i,j) + za(ip)*aaa
                                    z_alone_helper(i,j) = z_alone_helper(i,j) + aaa


                                    trans(a,b) = trans(a,b) +  aaa
                                    dummy(i,j) = dummy(i,j)  +  aaa
                                end if
                            else ! no offer made
                                call giveindex(i,j)
                                a=zx_idx !today index
                                call giveindex(ip,j)
                                b=zx_idx

                                mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*aaa
                                mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(ip,j)*aaa
                                mve_reg_helper(i,j) = mve_reg_helper(i,j) + aaa
                                z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*aaa
                                z_alone_end(i,j) = z_alone_end(i,j) + za(ip)*aaa
                                z_alone_helper(i,j) = z_alone_helper(i,j) + aaa



                                trans(a,b) = trans(a,b) +  aaa
                                dummy(i,j) = dummy(i,j)  +  aaa
                            end if
                        end do
                    end do
                end do
            end do


            ! meet -- perfect information case
            do ii = 1,numz ! other pre-innovaton productivity
                do jj = 1,numx ! other fixed cost

                    if (i>ii) then
                        prob_acq(i,j) = 1d0
                    elseif (i==ii) then
                        prob_acq(i,j) = 0.5
                    else
                        prob_acq(i,j) = 0d0
                    endif

                    do k = -1,1 ! my innovation realization
                        do m = -1,1 ! other innovation realization

                            ip = max(min(i+k,numz),1) ! target after innov
                            iip = max(min(ii+m,numz),1) ! acquirer after innov
                            inn_a = i_policy(ii,jj)
                            inn_t = i_policy(i,j)
                            aaa = Fs(ii,jj)/sum(Fs)*(1d0-info_asym)*(1d0-prob_acq(i,j)) *eta*p3(2-k,i,inn_t)*p3(2-m,ii,inn_a)*(1d0-psi)
                            if (sum(Fs) == 0d0) then
                                aaa = 0d0
                            end if


                            ! helper_var2 records the probability of a meeting conditional on target z and x

                            helper_var2(i,j) = helper_var2(i,j) + aaa

                            if (offermade_perfect(iip,jj,ip,j) == 1) then
                                ! helper_var3 records the probability of an offer being made conditional on target z and x and a meeting

                                helper_var3(i,j) = helper_var3(i,j) + aaa
                                helper_equity(i,j) = helper_equity(i,j) +  aaa*equity_ratio_perfect(iip,jj,ip,j)



                                if (indexs_perfect(iip,jj,ip,j)==1 .and. gc_perfect(iip,jj,ip,j)>1e-3) then
                                    new_corner_cash(i,j) = new_corner_cash(i,j) + aaa !;
                                end if

                                if (indexs_perfect(iip,jj,ip,j)/=1 .and. gc_perfect(iip,jj,ip,j)<=1e-3) then
                                    new_corner_equity(i,j) = new_corner_equity(i,j) + aaa !;
                                end if

                                if (indexs_perfect(iip,jj,ip,j)/=1 .and. gc_perfect(iip,jj,ip,j)>1e-3) then
                                    new_interior(i,j) = new_interior(i,j) + aaa !;
                                end if

                                if (indexs_perfect(iip,jj,ip,j)==1 .and. gc_perfect(iip,jj,ip,j)<=1e-3) then
                                    new_nonoffer(i,j) = new_nonoffer(i,j) + aaa!;
                                end if


                                if (gacp_eqm_perfect(iip,jj,ip,j) == 1) then ! if target accepts

                                    ! helper_var1 records the probability of acceptance conditional on target z and x, a meeting, and an offer made

                                    helper_var1(i,j) = helper_var1(i,j) + aaa
                                    helper_equity1(i,j) = helper_equity1(i,j) +  aaa*equity_ratio_perfect(iip,jj,ip,j)


                                    do ii_e = 1,numz ! entrant z
                                        do jj_e = 1,numx ! entrant x
                                            call giveindex(i,j)
                                            a=zx_idx !today index
                                            call giveindex(ii_e,jj_e)
                                            b=zx_idx

                                            F_exog_corrected = F_exog(ii_e,jj_e)/sum(F_exog)
                                            if (sum(F_exog)==0d0) then
                                                F_exog_corrected = 0d0
                                            end if


                                            trans(a,b) = trans(a,b) + aaa*(F_exog_corrected)
                                            dummy(i,j) = dummy(i,j) + aaa*(F_exog_corrected)
                                        end do
                                    end do
                                else ! if target rejects
                                    call giveindex(i,j)
                                    a=zx_idx !today index
                                    call giveindex(ip,j)
                                    b=zx_idx

                                    mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*aaa
                                    mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(ip,j)*aaa
                                    mve_reg_helper(i,j) = mve_reg_helper(i,j) + aaa
                                    z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*aaa
                                    z_alone_end(i,j) = z_alone_end(i,j) + za(ip)*aaa
                                    z_alone_helper(i,j) = z_alone_helper(i,j) + aaa


                                    trans(a,b) = trans(a,b) + aaa
                                    dummy(i,j) = dummy(i,j)  + aaa
                                end if
                            else ! no offer made
                                call giveindex(i,j)
                                a=zx_idx !today index
                                call giveindex(ip,j)
                                b=zx_idx

                                mve_reg_begin(i,j) = mve_reg_begin(i,j) + Vnew(i,j)*aaa
                                mve_reg_end(i,j) = mve_reg_end(i,j) + Vnew(ip,j)*aaa
                                mve_reg_helper(i,j) = mve_reg_helper(i,j) + aaa
                                z_alone_begin(i,j) = z_alone_begin(i,j) + za(i)*aaa
                                z_alone_end(i,j) = z_alone_end(i,j) + za(ip)*aaa
                                z_alone_helper(i,j) = z_alone_helper(i,j) + aaa


                                trans(a,b) = trans(a,b) + aaa
                                dummy(i,j) = dummy(i,j)  + aaa
                            end if
                        end do
                    end do
                end do
            end do
            !end if
        end do
    end do




    new_conditional_accept=helper_var1/helper_var3
    new_merger_prob=helper_var1 /helper_var2
    new_offer_prob=helper_var3/helper_var2
    new_equity_ratio_2=helper_equity/helper_var3
    new_equity_ratio_1=helper_equity1/helper_var3
    new_corner_cash = new_corner_cash/helper_var3
    new_corner_equity = new_corner_equity/helper_var3
    new_interior = new_interior/helper_var3
    new_nonoffer = new_nonoffer/helper_var3
    new_realized_gain = realized_gain/realized_gain_helper
    new_realized_gain_dollar = realized_gain_dollar/realized_gain_helper

    fixedcost_ratio=fixedcost_ratio/realized_gain_helper

    realized_gain_det = realized_gain_det/realized_gain_det_helper
    equity_ratio_acq = equity_ratio_acq/realized_gain_det_helper
    value_loss = value_loss/realized_gain_det_helper
    value_loss2 = value_loss2/realized_gain_det_helper

    tar_frac_of_gain = tar_frac_of_gain/realized_gain_helper

    acq_ann_ret = acq_ann_ret/acq_ann_ret_helper
    tar_ann_ret = tar_ann_ret/acq_ann_ret_helper
    tar_frac_merger_gain = tar_frac_merger_gain/acq_ann_ret_helper
    cash_ratio =  cash_ratio/acq_ann_ret_helper
    log_mv_ratio = log_mv_ratio/acq_ann_ret_helper
    rel_size = rel_size/acq_ann_ret_helper

    mve_reg_end = mve_reg_end/mve_reg_helper
    mve_reg_begin = mve_reg_begin/mve_reg_helper

    z_alone_begin = z_alone_begin/z_alone_helper
    z_alone_end = z_alone_end/z_alone_helper


    avg_rnd_intensity = 0d0
    avg_innovation = 0d0
    avg_i_policy = 0d0
    avg_growth_rate = 0d0
    agg_rnd=0d0
    do i=1,numz
        do j=1,numx
            avg_rnd_intensity = avg_rnd_intensity + (Fs(i,j)/sum(Fs))*(i_cost(i_policy(i,j))*(zeta/pi))
            avg_innovation = avg_innovation + (Fs(i,j)/sum(Fs))*(mu_value(i_policy(i,j)))
            avg_i_policy = avg_i_policy + (Fs(i,j)/sum(Fs))*(real(i_policy(i,j)))
            avg_growth_rate = avg_growth_rate + (Fs(i,j)/sum(Fs))*(z_alone_end(i,j)/(z_alone_end(i,j)+z_alone_begin(i,j))*2d0-1d0)
            agg_rnd= agg_rnd+ (Fs(i,j)/sum(Fs))*(i_cost(i_policy(i,j)))*za(i)
        enddo
    enddo

    growth_store=z_alone_end/(z_alone_end+z_alone_begin)*2d0-1d0



    do i=1,numz
        do j=1,numx
            Ft(i,j) = Fs(i,j) * helper_var3(i,j)
        enddo
    enddo

    ! Normalizations
    mean_z = 0d0
    do i=1,numz
        do j=1,numx
            mean_z = mean_z + (Fs(i,j)/sum(Fs))*log(zt(i))
        enddo
    enddo
    dummy_z = 0d0
    do i=1,numz
        do j=1,numx
            dummy_z = dummy_z + (Fs(i,j)/sum(Fs))*((log(zt(i)) - mean_z)**2d0)
        enddo
    enddo
    std_z = sqrt(dummy_z)
    normalized_z = ((log(zt) - mean_z)/std_z)* 2.098416 + 4.792275





    mean_zt = 0d0
    do i=1,numz
        do j=1,numx
            mean_zt = mean_zt + (Ft(i,j)/sum(Ft))*log(zt(i))
        enddo
    enddo
    dummy_zt = 0d0
    do i=1,numz
        do j=1,numx
            dummy_zt = dummy_zt + (Ft(i,j)/sum(Ft))*((log(zt(i)) - mean_zt)**2d0)
        enddo
    enddo
    std_zt = sqrt(dummy_zt)
    normalized_zt = ((log(zt) - mean_zt)/std_zt)*  1.924635 + 5.028083   ! For percentage of stock regression
    normalized_zt2 = ((log(zt) - mean_zt)/std_zt)*  1.901221 +  4.995369   ! For deal completion regression



    avg_value_loss = 0d0
    avg_value_loss_helper = 0d0
    agg_value_loss= 0d0
    do ip=1,numz
        do j=1,numx
            do iip=1,numz
                do jj=1,numx
                    do i=1,2 ! interp realization
                        do ii=1,2 ! kesi realization
                            if (isnan(value_loss(ip,j,iip,jj,i,ii))==.FALSE.) then
                                avg_value_loss = avg_value_loss + value_loss(ip,j,iip,jj,i,ii)*realized_gain_det_helper(ip,j,iip,jj,i,ii)
                                agg_value_loss = agg_value_loss + value_loss2(ip,j,iip,jj,i,ii)*realized_gain_det_helper(ip,j,iip,jj,i,ii)
                                avg_value_loss_helper = avg_value_loss_helper + realized_gain_det_helper(ip,j,iip,jj,i,ii)
                            endif
                        enddo
                    enddo
                enddo
            enddo
        enddo
    enddo
    avg_value_loss = avg_value_loss/avg_value_loss_helper
    agg_value_loss = agg_value_loss/avg_value_loss_helper




    avg_acq_ann_ret = 0d0
    avg_acq_ann_ret_helper = 0d0
    do i=1,numz
        do j=1,numx
            do ii=1,numz
                do jj=1,numx

                    if (isnan(acq_ann_ret(i,j,ii,jj))==.FALSE.) then
                        avg_acq_ann_ret = avg_acq_ann_ret + acq_ann_ret(i,j,ii,jj)*acq_ann_ret_helper(i,j,ii,jj)*Fs(i,j)
                        avg_acq_ann_ret_helper = avg_acq_ann_ret_helper + acq_ann_ret_helper(i,j,ii,jj)*Fs(i,j)
                    endif

                enddo
            enddo
        enddo
    enddo
    avg_acq_ann_ret = avg_acq_ann_ret/avg_acq_ann_ret_helper

    avg_tar_ann_ret = 0d0
    avg_tar_ann_ret_helper = 0d0
    do i=1,numz
        do j=1,numx
            do ii=1,numz
                do jj=1,numx

                    if (isnan(acq_ann_ret(i,j,ii,jj))==.FALSE.) then
                        avg_tar_ann_ret = avg_tar_ann_ret + tar_ann_ret(i,j,ii,jj)*acq_ann_ret_helper(i,j,ii,jj)*Fs(i,j)
                        avg_tar_ann_ret_helper = avg_tar_ann_ret_helper + acq_ann_ret_helper(i,j,ii,jj)*Fs(i,j)
                    endif

                enddo
            enddo
        enddo
    enddo
    avg_tar_ann_ret = avg_tar_ann_ret/avg_tar_ann_ret_helper



    avg_tar_frac_merger_gain = 0d0
    avg_tar_frac_merger_gain_helper = 0d0
    do i=1,numz
        do j=1,numx
            do ii=1,numz
                do jj=1,numx

                    if (isnan(acq_ann_ret(i,j,ii,jj))==.FALSE.) then
                        avg_tar_frac_merger_gain = avg_tar_frac_merger_gain + tar_frac_merger_gain(i,j,ii,jj)*acq_ann_ret_helper(i,j,ii,jj)*Fs(i,j)
                        avg_tar_frac_merger_gain_helper = avg_tar_frac_merger_gain_helper + acq_ann_ret_helper(i,j,ii,jj)*Fs(i,j)
                    endif

                enddo
            enddo
        enddo
    enddo
    avg_tar_frac_merger_gain = avg_tar_frac_merger_gain/avg_tar_frac_merger_gain_helper

    avg_rel_size = 0d0
    avg_rel_size_helper = 0d0
    do i=1,numz
        do j=1,numx
            do ii=1,numz
                do jj=1,numx

                    if (isnan(acq_ann_ret(i,j,ii,jj))==.FALSE..and. rel_size(i,j,ii,jj)/=0) then
                        avg_rel_size = avg_rel_size + rel_size(i,j,ii,jj)*acq_ann_ret_helper(i,j,ii,jj)*Fs(i,j)
                        avg_rel_size_helper = avg_rel_size_helper + acq_ann_ret_helper(i,j,ii,jj)*Fs(i,j)
                    endif

                enddo
            enddo
        enddo
    enddo
    avg_rel_size = avg_rel_size/avg_rel_size_helper
  



    F_acq=Fs*realized_gain_helper
    F_acq=F_acq/sum(F_acq)
    F_tar=Fs*helper_var1
    F_tar=F_tar/sum(F_tar)





    avg_merger_prob = 0d0
    do i=1,numz
        do j=1,numx
            if (isnan(new_merger_prob(i,j))==.FALSE.) then
                avg_merger_prob = avg_merger_prob + 2d0*helper_var1(i,j)*Fs(i,j)
            endif
        enddo
    enddo



    avg_realized_gain = 0d0
    avg_fixedcost_ratio = 0d0
    avg_tar_frac_of_gain=0d0
    avg_realized_gain_helper = 0d0
    avg_realized_gain_dollar=0d0
    avg_realized_gain_dollar_unconditional=0d0
    do i=1,numz
        do j=1,numx
            if (isnan(new_realized_gain(i,j))==.FALSE.) then
                avg_realized_gain = avg_realized_gain + new_realized_gain(i,j)*Fs(i,j)   * realized_gain_helper(i,j)
                avg_realized_gain_dollar = avg_realized_gain_dollar + new_realized_gain_dollar(i,j)*Fs(i,j)   * realized_gain_helper(i,j)

                avg_realized_gain_helper = avg_realized_gain_helper + Fs(i,j) * realized_gain_helper(i,j)

                avg_fixedcost_ratio = avg_fixedcost_ratio + fixedcost_ratio(i,j)*Fs(i,j) * realized_gain_helper(i,j)

                avg_tar_frac_of_gain = avg_tar_frac_of_gain + tar_frac_of_gain(i,j)*Fs(i,j)*realized_gain_helper(i,j)
            endif
        enddo
    enddo
    avg_realized_gain_dollar_unconditional = avg_realized_gain_dollar
    avg_realized_gain = avg_realized_gain/avg_realized_gain_helper
    avg_realized_gain_dollar = avg_realized_gain_dollar/avg_realized_gain_helper







    avg_fixedcost_ratio = avg_fixedcost_ratio/avg_realized_gain_helper
    avg_tar_frac_of_gain = avg_tar_frac_of_gain/avg_realized_gain_helper






    e_ratio = 0d0
    zt_value = 0d0
    x_value = 0d0
    omega_value = 0d0

    ii=1
    do i=1,numz
        do j=1,numx
            if (isnan(new_equity_ratio_2(i,j))==.FALSE.) then

                e_ratio(ii,1)=new_equity_ratio_2(i,j)
                zt_value(ii,1)=normalized_zt(i)
                !x_value(ii,1)=x(j)
                omega_value(ii,1)=Ft(i,j)
                ii=ii+1
            end if
        enddo
    enddo

    allocate(Xmatrix(ii-1,2))
    allocate(Xmatrixt(2,ii-1))
    allocate(XWX(2,2))
    allocate(XWXinv(2,2))
    allocate(omegaw(ii-1,ii-1))
    allocate(Ymatrix(ii-1,1))
    allocate(XWY(2,1))
    allocate(y_vector(ii-1,1))

    omegaw=0
    do jj=1,ii-1
        omegaw(jj,jj) = omega_value(jj,1)
    enddo
    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=zt_value(1:ii-1,1)
    y_vector(:,1)=e_ratio(1:ii-1,1)

    Xmatrixt=transpose(Xmatrix)
    Xmatrix=matmul(omegaw,Xmatrix)
    XWX=matmul(Xmatrixt,Xmatrix)

    call findinv(XWX, XWXinv, 2, errorflag)

    Ymatrix=matmul(omegaw,y_vector)
    XWY=matmul(Xmatrixt,Ymatrix)

    bbeta=matmul(XWXinv,XWY)


    mean_equity_ratio=sum(e_ratio*omega_value/sum(omega_value))



    deallocate(Xmatrix)
    deallocate(Xmatrixt)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(omegaw)
    deallocate(Ymatrix)
    deallocate(XWY)
    deallocate(y_vector)



    e_ratio = 0d0
    zt_value = 0d0
    x_value = 0d0
    omega_value = 0d0

    new_merger_prob=helper_var1
    new_merger_prob_unconditional=new_merger_prob
    ii=1
    do i=1,numz
        do j=1,numx
            if (isnan(new_merger_prob(i,j))==.FALSE.) then

                e_ratio(ii,1)=new_merger_prob(i,j)
                zt_value(ii,1)=normalized_z(i)
                omega_value(ii,1)=Fs(i,j)
                ii=ii+1
            end if
        enddo
    enddo

    allocate(Xmatrix(ii-1,3))
    allocate(Xmatrixt(3,ii-1))
    allocate(XWX(3,3))
    allocate(XWXinv(3,3))
    allocate(omegaw(ii-1,ii-1))
    allocate(Ymatrix(ii-1,1))
    allocate(XWY(3,1))
    allocate(y_vector(ii-1,1))

    omegaw=0
    do jj=1,ii-1
        omegaw(jj,jj) = omega_value(jj,1)
    enddo

    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=zt_value(1:ii-1,1)
    Xmatrix(:,3)=zt_value(1:ii-1,1)**2d0
    y_vector(:,1)=e_ratio(1:ii-1,1)



    Xmatrixt=transpose(Xmatrix)
    Xmatrix=matmul(omegaw,Xmatrix)
    XWX=matmul(Xmatrixt,Xmatrix)

    call findinv(XWX, XWXinv, 3, errorflag)

    Ymatrix=matmul(omegaw,y_vector)
    XWY=matmul(Xmatrixt,Ymatrix)

    bbeta2=matmul(XWXinv,XWY)

    !!!!!!!

    deallocate(Xmatrix)
    deallocate(Xmatrixt)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(omegaw)
    deallocate(Ymatrix)
    deallocate(XWY)
    deallocate(y_vector)

    e_ratio = 0d0
    zt_value = 0d0
    x_value = 0d0
    omega_value = 0d0

    new_merger_prob=helper_var1/helper_var2

    ii=1
    do i=1,numz
        do j=1,numx
            if (isnan(new_merger_prob(i,j))==.FALSE.) then

                e_ratio(ii,1)=new_merger_prob(i,j)
                zt_value(ii,1)=normalized_zt(i)
                !x_value(ii,1)=x(j)
                omega_value(ii,1)=Ft(i,j)
                ii=ii+1
            end if
        enddo
    enddo

    allocate(Xmatrix(ii-1,3))
    allocate(Xmatrixt(3,ii-1))
    allocate(XWX(3,3))
    allocate(XWXinv(3,3))
    allocate(omegaw(ii-1,ii-1))
    allocate(Ymatrix(ii-1,1))
    allocate(XWY(3,1))
    allocate(y_vector(ii-1,1))

    omegaw=0
    do jj=1,ii-1
        omegaw(jj,jj) = omega_value(jj,1)
    enddo

    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=zt_value(1:ii-1,1)
    Xmatrix(:,3)=zt_value(1:ii-1,1)**2d0
    y_vector(:,1)=e_ratio(1:ii-1,1)

    Xmatrixt=transpose(Xmatrix)
    Xmatrix=matmul(omegaw,Xmatrix)
    XWX=matmul(Xmatrixt,Xmatrix)

    call findinv(XWX, XWXinv, 3, errorflag)

    Ymatrix=matmul(omegaw,y_vector)
    XWY=matmul(Xmatrixt,Ymatrix)

    bbeta4=matmul(XWXinv,XWY)



    zbar = 0d0
    do i=1,numz
        do j=1,numx
            zbar = zbar + zt(i)*Fs(i,j)
        enddo
    enddo

    output = pi/zeta*zbar
    agg_rnd_intensity=agg_rnd/output


    consumption=output-agg_rnd-(agg_value_loss*avg_merger_prob/2d0) -delta*kappa/(r+delta)*output



    deallocate(Xmatrix)
    deallocate(Xmatrixt)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(omegaw)
    deallocate(Ymatrix)
    deallocate(XWY)
    deallocate(y_vector)

    e_ratio = 0d0
    zt_value = 0d0
    x_value = 0d0
    omega_value = 0d0

    ii=1
    do i=1,numz
        do j=1,numx
            if (isnan(new_conditional_accept(i,j))==.FALSE.) then

                e_ratio(ii,1)=new_conditional_accept(i,j)
                zt_value(ii,1)=normalized_zt2(i)
                !x_value(ii,1)=x(j)
                omega_value(ii,1)=Ft(i,j)
                ii=ii+1
            end if
        enddo
    enddo

    allocate(Xmatrix(ii-1,2))
    allocate(Xmatrixt(2,ii-1))
    allocate(XWX(2,2))
    allocate(XWXinv(2,2))
    allocate(omegaw(ii-1,ii-1))
    allocate(Ymatrix(ii-1,1))
    allocate(XWY(2,1))
    allocate(y_vector(ii-1,1))

    omegaw=0
    do jj=1,ii-1
        omegaw(jj,jj) = omega_value(jj,1)
    enddo

    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=zt_value(1:ii-1,1)
    y_vector(:,1)=e_ratio(1:ii-1,1)

    Xmatrixt=transpose(Xmatrix)
    Xmatrix=matmul(omegaw,Xmatrix)
    XWX=matmul(Xmatrixt,Xmatrix)

    call findinv(XWX, XWXinv, 2, errorflag)

    Ymatrix=matmul(omegaw,y_vector)
    XWY=matmul(Xmatrixt,Ymatrix)

    bbeta3=matmul(XWXinv,XWY)



    deallocate(Xmatrix)
    deallocate(Xmatrixt)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(omegaw)
    deallocate(Ymatrix)
    deallocate(XWY)
    deallocate(y_vector)



    e_ratio = 0d0
    zt_value = 0d0
    x_value = 0d0
    omega_value = 0d0

    ii=1
    do i=1,numz
        do j=1,numx
            do mm=1,numz
                do jj=1,numx
                    if (isnan(acq_ann_ret(i,j,mm,jj))==.FALSE. ) then
                        e_ratio(ii,1)=acq_ann_ret(i,j,mm,jj)
                        zt_value(ii,1)=cash_ratio(i,j,mm,jj)
                        x_value(ii,1)=log_mv_ratio(i,j,mm,jj)
                        omega_value(ii,1)=acq_ann_ret_helper(i,j,mm,jj)*Fs(i,j) ! *Fs(mm,jj)
                        zt_value2(ii,1)=normalized_zt2(mm)
                        ii=ii+1
                    end if
                enddo
            enddo
        enddo
    enddo



    allocate(Xmatrix(ii-1,4))
    allocate(Xmatrixt(4,ii-1))
    allocate(XWX(4,4))
    allocate(XWXinv(4,4))
    allocate(omegaw(ii-1,ii-1))
    allocate(Ymatrix(ii-1,1))
    allocate(XWY(4,1))
    allocate(y_vector(ii-1,1))

    omegaw=0
    do jj=1,ii-1
        omegaw(jj,jj) = omega_value(jj,1)
    enddo

    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=zt_value(1:ii-1,1)
    Xmatrix(:,3)=x_value(1:ii-1,1)
    Xmatrix(:,4)=zt_value2(1:ii-1,1)
    y_vector(:,1)=e_ratio(1:ii-1,1)



    Xmatrixt=transpose(Xmatrix)
    Xmatrix=matmul(omegaw,Xmatrix)
    XWX=matmul(Xmatrixt,Xmatrix)

    call findinv(XWX, XWXinv, 4, errorflag)

    Ymatrix=matmul(omegaw,y_vector)
    XWY=matmul(Xmatrixt,Ymatrix)

    bbeta5=matmul(XWXinv,XWY)

    deallocate(Xmatrix)
    deallocate(Xmatrixt)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(omegaw)
    deallocate(Ymatrix)
    deallocate(XWY)
    deallocate(y_vector)

    e_ratio = 0d0
    zt_value = 0d0
    x_value = 0d0
    omega_value = 0d0

    ii=1
    do i=1,numz
        do j=1,numx
            if (isnan(realized_gain_helper(i,j))==.FALSE.) then

                e_ratio(ii,1)=realized_gain_helper(i,j)
                zt_value(ii,1)=normalized_z(i)
                !x_value(ii,1)=x(j)
                omega_value(ii,1)=Fs(i,j)
                ii=ii+1
            end if
        enddo
    enddo

    allocate(Xmatrix(ii-1,2))
    allocate(Xmatrixt(2,ii-1))
    allocate(XWX(2,2))
    allocate(XWXinv(2,2))
    allocate(omegaw(ii-1,ii-1))
    allocate(Ymatrix(ii-1,1))
    allocate(XWY(2,1))
    allocate(y_vector(ii-1,1))

    omegaw=0
    do jj=1,ii-1
        omegaw(jj,jj) = omega_value(jj,1)
    enddo
    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=zt_value(1:ii-1,1)
    y_vector(:,1)=e_ratio(1:ii-1,1)

    Xmatrixt=transpose(Xmatrix)
    Xmatrix=matmul(omegaw,Xmatrix)
    XWX=matmul(Xmatrixt,Xmatrix)

    call findinv(XWX, XWXinv, 2, errorflag)

    Ymatrix=matmul(omegaw,y_vector)
    XWY=matmul(Xmatrixt,Ymatrix)

    bbeta6=matmul(XWXinv,XWY)

    deallocate(Xmatrix)
    deallocate(Xmatrixt)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(omegaw)
    deallocate(Ymatrix)
    deallocate(XWY)
    deallocate(y_vector)


    e_ratio = 0d0
    zt_value = 0d0
    x_value = 0d0
    omega_value = 0d0

    ii=1
    do i=1,numz
        do j=1,numx
            if (isnan(mve_reg_end(i,j))==.FALSE.) then
                e_ratio(ii,1)=log(mve_reg_end(i,j))
                zt_value(ii,1)=log(mve_reg_begin(i,j))
                omega_value(ii,1)=Fs(i,j)*mve_reg_helper(i,j)
                ii=ii+1
            end if
        enddo
    enddo

    allocate(Xmatrix(ii-1,2))
    allocate(Xmatrixt(2,ii-1))
    allocate(XWX(2,2))
    allocate(XWXinv(2,2))
    allocate(omegaw(ii-1,ii-1))
    allocate(Ymatrix(ii-1,1))
    allocate(XWY(2,1))
    allocate(y_vector(ii-1,1))

    allocate(eps_res(ii-1,1))

    omegaw=0
    do jj=1,ii-1
        omegaw(jj,jj) = omega_value(jj,1)
    enddo

    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=zt_value(1:ii-1,1)
    mean_log_mve = sum(zt_value(1:ii-1,1)*omega_value(1:ii-1,1))/sum(omega_value(1:ii-1,1))
    std_log_mve = sqrt((sum(((zt_value(1:ii-1,1)-mean_log_mve)**2d0)*omega_value(1:ii-1,1)))/sum(omega_value(1:ii-1,1)))
   
    Xmatrix(:,2) = (Xmatrix(:,2) - mean_log_mve)  / std_log_mve

    y_vector(:,1)=e_ratio(1:ii-1,1)
    y_vector(:,1)= (y_vector(:,1)-mean_log_mve)  / std_log_mve

    Xmatrixt=transpose(Xmatrix)
    Xmatrix=matmul(omegaw,Xmatrix)
    XWX=matmul(Xmatrixt,Xmatrix)

    call findinv(XWX, XWXinv, 2, errorflag)

    Ymatrix=matmul(omegaw,y_vector)
    XWY=matmul(Xmatrixt,Ymatrix)

    bbeta7=matmul(XWXinv,XWY)


    eps_res(:,1) =  y_vector(:,1)-matmul(transpose(Xmatrixt),bbeta7(:,1))
    std_eps_res = sqrt(sum(((eps_res(1:ii-1,1))**2d0)*omega_value(1:ii-1,1))/sum(omega_value(1:ii-1,1)))


    deallocate(Xmatrix)
    deallocate(Xmatrixt)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(omegaw)
    deallocate(Ymatrix)
    deallocate(XWY)
    deallocate(y_vector)
    deallocate(eps_res)

    e_ratio = 0d0
    zt_value = 0d0
    x_value = 0d0
    omega_value = 0d0

    ii=1
    do i=1,numz
        do j=1,numx
            do mm=1,numz
                do jj=1,numx
                    if (isnan(acq_ann_ret(i,j,mm,jj))==.FALSE. ) then
                        e_ratio(ii,1)=acq_ann_ret(i,j,mm,jj)
                        zt_value(ii,1)=cash_ratio(i,j,mm,jj)
                        x_value(ii,1)=log_mv_ratio(i,j,mm,jj)
                        omega_value(ii,1)=acq_ann_ret_helper(i,j,mm,jj)*Fs(i,j) ! *Fs(mm,jj)
                        ii=ii+1
                    end if
                enddo
            enddo
        enddo
    enddo

    allocate(Xmatrix(ii-1,3))
    allocate(Xmatrixt(3,ii-1))
    allocate(XWX(3,3))
    allocate(XWXinv(3,3))
    allocate(omegaw(ii-1,ii-1))
    allocate(Ymatrix(ii-1,1))
    allocate(XWY(3,1))
    allocate(y_vector(ii-1,1))



    omegaw=0
    do jj=1,ii-1
        omegaw(jj,jj) = omega_value(jj,1)
    enddo

    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=zt_value(1:ii-1,1)
    Xmatrix(:,3)=x_value(1:ii-1,1)
    y_vector(:,1)=e_ratio(1:ii-1,1)





    Xmatrixt=transpose(Xmatrix)
    Xmatrix=matmul(omegaw,Xmatrix)
    XWX=matmul(Xmatrixt,Xmatrix)

    call findinv(XWX, XWXinv, 3, errorflag)

    Ymatrix=matmul(omegaw,y_vector)
    XWY=matmul(Xmatrixt,Ymatrix)

    bbeta8=matmul(XWXinv,XWY)

    deallocate(Xmatrix)
    deallocate(Xmatrixt)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(omegaw)
    deallocate(Ymatrix)
    deallocate(XWY)
    deallocate(y_vector)








    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Simulation!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call cpu_time(start)


    nsim1=20000
    nt1=27 !27

    allocate(initial_z_shock(nsim1,1))
    allocate(psi_shock(nsim1,nt1))
    allocate(eta_shock(nsim1,nt1))
    allocate(innov_shock(nsim1,nt1))
    allocate(innov_other_shock(nsim1,nt1))
    allocate(z_other_shock(nsim1,nt1))
    allocate(tiebreaking_shock(nsim1,nt1))
    allocate(omega_shock(nsim1,nt1))
    allocate(merger_interp_shock(nsim1,nt1))




    m=nsim1*nt1*10

    call random_seed(SIZE = m)      
    call random_number(initial_z_shock)
    call random_number(psi_shock)
    call random_number(eta_shock)
    call random_number(innov_shock)
    call random_number(innov_other_shock)
    call random_number(z_other_shock)
    call random_number(tiebreaking_shock)
    call random_number(omega_shock)
    call random_number(merger_interp_shock)


    if (estimate==0.and.comparastat==0) then

        OPEN(1,file='initial_z_shock.txt') ; WRITE(1,"(f20.8)") initial_z_shock
        OPEN(1,file='psi_shock.txt') ; WRITE(1,"(f20.8)") psi_shock
        OPEN(1,file='eta_shock.txt') ; WRITE(1,"(f20.8)") eta_shock
        OPEN(1,file='innov_shock.txt') ; WRITE(1,"(f20.8)") innov_shock
        OPEN(1,file='innov_other_shock.txt') ; WRITE(1,"(f20.8)") innov_other_shock
        OPEN(1,file='z_other_shock.txt') ; WRITE(1,"(f20.8)") z_other_shock
        OPEN(1,file='tiebreaking_shock.txt') ; WRITE(1,"(f20.8)") tiebreaking_shock
        OPEN(1,file='omega_shock.txt') ; WRITE(1,"(f20.8)") omega_shock
        OPEN(1,file='merger_interp_shock.txt') ; WRITE(1,"(f20.8)") merger_interp_shock
    end if

    !print *, 'initial_z_shock(1:10)',initial_z_shock(1:10,1)
    !print *, 'eta_shock(1:10)',eta_shock(1:10,1)
    !print *, 'omega_shock(1:10)',omega_shock(1:10,1)
    !pause

    allocate(sim_owner(nsim1,nt1+1))
    allocate(sim_z(nsim1,nt1+1))
    allocate(sim_z_preinnov(nsim1,nt1))
    allocate(sim_z_postinnov(nsim1,nt1))
    allocate(sim_was_acquired_dummy(nsim1,nt1))
    allocate(sim_rec_offer_dummy(nsim1,nt1))
    allocate(sim_has_acquired_dummy(nsim1,nt1))
    allocate(sim_equity_ratio(nsim1,nt1))
   

    Fscum=0d0
    do i = 2,numz+1
        Fscum(i) = sum(Fs(1:i-1,1))
    enddo


    sim_owner=0
    sim_z=0
    sim_z_preinnov=0
    sim_z_postinnov=0
    sim_was_acquired_dummy=0
    sim_rec_offer_dummy=0
    sim_has_acquired_dummy=0
    sim_equity_ratio=-1d0

    do n=1,nsim1
        do i=1,numz
            if ((initial_z_shock(n,1)>=Fscum(i)).and.(initial_z_shock(n,1)<Fscum(i+1))) then
                sim_z(n,1)=i
            endif
        enddo
    enddo

    sim_owner(:,1) = 1


    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SIMULATION CODE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    p3_cum = 0d0 
    do j = 1,3 
        do i = 1,numz
            do inn = 1,numi
                p3_cum(j,i,inn) = sum(p3(1:j,i,inn))
            enddo
        enddo
    enddo


    j = 1
    jj = 1



    do n = 1,nsim1
        do t  = 1,nt1
            if (psi_shock(n,t)<psi) then ! exogenous exit

                sim_owner(n,t+1)=sim_owner(n,t)+1
                sim_z(n,t+1)=1
                sim_z_preinnov(n,t)= -1
                sim_z_postinnov(n,t)= -1
                sim_was_acquired_dummy(n,t)= -1
                sim_rec_offer_dummy(n,t)= -1
                sim_has_acquired_dummy(n,t)= -1

            else ! no exogenous exit

                i = sim_z(n,t)

                if (eta_shock(n,t)>eta) then ! no meeting

                    if (innov_shock(n,t) < p3_cum(1,i,i_policy(i,j))) then ! high z'

                        sim_owner(n,t+1)=sim_owner(n,t)
                        sim_z(n,t+1)= min(i+1,numz)
                        sim_z_preinnov(n,t)= i
                        sim_z_postinnov(n,t)= min(i+1,numz)
                        sim_was_acquired_dummy(n,t)= 0
                        sim_rec_offer_dummy(n,t)= 0
                        sim_has_acquired_dummy(n,t)= 0

                    elseif (innov_shock(n,t) < p3_cum(2,i,i_policy(i,j))) then ! med z'

                        sim_owner(n,t+1)=sim_owner(n,t)
                        sim_z(n,t+1)= i
                        sim_z_preinnov(n,t)= i
                        sim_z_postinnov(n,t)= i
                        sim_was_acquired_dummy(n,t)= 0
                        sim_rec_offer_dummy(n,t)= 0
                        sim_has_acquired_dummy(n,t)= 0

                    else ! low z'

                        sim_owner(n,t+1)=sim_owner(n,t)
                        sim_z(n,t+1)= max(i-1,1)
                        sim_z_preinnov(n,t)= i
                        sim_z_postinnov(n,t)= max(i-1,1)
                        sim_was_acquired_dummy(n,t)= 0
                        sim_rec_offer_dummy(n,t)= 0
                        sim_has_acquired_dummy(n,t)= 0

                    endif

                else ! meeting

                    do ii_draw=1,numz ! draw the z of the met firm
                        if ((z_other_shock(n,t)>=Fscum(ii_draw)).and.(z_other_shock(n,t)<Fscum(ii_draw+1))) then
                            ii=ii_draw
                        endif
                    enddo

                    if (i>ii) then ! determine who is the acquirer and who is the target
                        is_acquirer = 1
                    elseif (i==ii) then
                        if (tiebreaking_shock(n,t) < 0.5) then
                            is_acquirer = 1
                        else
                            is_acquirer = 0
                        endif
                    else
                        is_acquirer = 0
                    endif

                    if (omega_shock(n,t) < info_asym) then
                        imperfect_info = 1
                    else
                        imperfect_info = 0
                    endif

                    ! determine the new productivity of our firm
                    if (innov_shock(n,t) < p3_cum(1,i,i_policy(i,j))) then ! high z'
                        ip = min(i+1,numz)
                        k = 1
                    elseif (innov_shock(n,t) < p3_cum(2,i,i_policy(i,j))) then ! med z'
                        ip = i
                        k = 0
                    else ! low z'
                        ip = max(i-1,1)
                        k = -1
                    endif

                    ! determine the new productivity of the other firm
                    if (innov_shock(n,t) < p3_cum(1,ii,i_policy(ii,jj))) then ! high z'
                        iip = min(ii+1,numz)
                        m = 1
                    elseif (innov_shock(n,t) < p3_cum(2,ii,i_policy(ii,jj))) then ! med z'
                        iip = ii
                        m = 0
                    else ! low z'
                        iip = max(ii-1,1)
                        m = -1
                    endif

                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    if ((imperfect_info == 1) .and. (is_acquirer == 1)) then ! is acquirer and imperfect info
                        if (offermade_i(ip,j,ii,jj,i_policy(ii,jj)) == 1) then ! if we make an offer
                            if (gacp_eqm_i(ip,j,ii,jj,2-m,i_policy(ii,jj)) == 1) then ! if the target accepts

                                ![p_lb,p_ub,i_lb,i_ub]=zinterp(ip,iip)
                                call zinterp(ip,iip,p_lb, p_ub,i_lb,i_ub)

                                if (merger_interp_shock(n,t)<p_lb) then ! lower case is realized

                                    sim_owner(n,t+1)=sim_owner(n,t)
                                    sim_z(n,t+1)= i_lb
                                    sim_z_preinnov(n,t)= i
                                    sim_z_postinnov(n,t)= ip
                                    sim_was_acquired_dummy(n,t)= 0
                                    sim_rec_offer_dummy(n,t)= 0
                                    sim_has_acquired_dummy(n,t)= 1

                                else ! upper case is realized

                                    sim_owner(n,t+1)=sim_owner(n,t)
                                    sim_z(n,t+1)= i_ub
                                    sim_z_preinnov(n,t)= i
                                    sim_z_postinnov(n,t)= ip
                                    sim_was_acquired_dummy(n,t)= 0
                                    sim_rec_offer_dummy(n,t)= 0
                                    sim_has_acquired_dummy(n,t)= 1

                                endif

                            else ! if we make an offer, but the target rejects

                                sim_owner(n,t+1)=sim_owner(n,t)
                                sim_z(n,t+1)= ip
                                sim_z_preinnov(n,t)= i
                                sim_z_postinnov(n,t)= ip
                                sim_was_acquired_dummy(n,t)= 0
                                sim_rec_offer_dummy(n,t)= 0
                                sim_has_acquired_dummy(n,t)= 0

                            endif
                        else ! if we make no offer

                            sim_owner(n,t+1)=sim_owner(n,t)
                            sim_z(n,t+1)= ip
                            sim_z_preinnov(n,t)= i
                            sim_z_postinnov(n,t)= ip
                            sim_was_acquired_dummy(n,t)= 0
                            sim_rec_offer_dummy(n,t)= 0
                            sim_has_acquired_dummy(n,t)= 0

                        endif
                        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    elseif ((imperfect_info == 0) .and. (is_acquirer == 1)) then ! is acquirer and perfect info
                        if (offermade_perfect(ip,j,iip,jj) == 1) then ! if we make an offer
                            if (gacp_eqm_perfect(ip,j,iip,jj) == 1) then ! if the target accepts

                                ![p_lb,p_ub,i_lb,i_ub]=zinterp(ip,iip)
                                call zinterp(ip,iip,p_lb, p_ub,i_lb,i_ub)

                                if (merger_interp_shock(n,t)<p_lb) then ! lower case is realized

                                    sim_owner(n,t+1)=sim_owner(n,t)
                                    sim_z(n,t+1)= i_lb
                                    sim_z_preinnov(n,t)= i
                                    sim_z_postinnov(n,t)= ip
                                    sim_was_acquired_dummy(n,t)= 0
                                    sim_rec_offer_dummy(n,t)= 0
                                    sim_has_acquired_dummy(n,t)= 1

                                else ! upper case is realized

                                    sim_owner(n,t+1)=sim_owner(n,t)
                                    sim_z(n,t+1)= i_ub
                                    sim_z_preinnov(n,t)= i
                                    sim_z_postinnov(n,t)= ip
                                    sim_was_acquired_dummy(n,t)= 0
                                    sim_rec_offer_dummy(n,t)= 0
                                    sim_has_acquired_dummy(n,t)= 1

                                endif

                            else ! if we make an offer, but the target rejects

                                sim_owner(n,t+1)=sim_owner(n,t)
                                sim_z(n,t+1)= ip
                                sim_z_preinnov(n,t)= i
                                sim_z_postinnov(n,t)= ip
                                sim_was_acquired_dummy(n,t)= 0
                                sim_rec_offer_dummy(n,t)= 0
                                sim_has_acquired_dummy(n,t)= 0

                            endif
                        else ! if we make no offer

                            sim_owner(n,t+1)=sim_owner(n,t)
                            sim_z(n,t+1)= ip
                            sim_z_preinnov(n,t)= i
                            sim_z_postinnov(n,t)= ip
                            sim_was_acquired_dummy(n,t)= 0
                            sim_rec_offer_dummy(n,t)= 0
                            sim_has_acquired_dummy(n,t)= 0

                        endif
                        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    elseif ((imperfect_info == 1) .and. (is_acquirer == 0)) then ! is target and imperfect info
                        if (offermade_i(iip,jj,i,j,i_policy(i,j)) == 1) then ! if the other firm made us an offer
                            if (gacp_eqm_i(iip,jj,i,j,2-k,i_policy(i,j)) == 1) then  ! if we accept the offer

                                sim_owner(n,t+1)=sim_owner(n,t)+1
                                sim_z(n,t+1)=1
                                sim_z_preinnov(n,t)= i
                                sim_z_postinnov(n,t)= ip
                                sim_was_acquired_dummy(n,t)= 1
                                sim_rec_offer_dummy(n,t)= 1
                                sim_has_acquired_dummy(n,t)= 0
                                sim_equity_ratio(n,t)=equity_ratio(iip,jj,i,j)

                            else ! if we reject the offer

                                sim_owner(n,t+1)=sim_owner(n,t)
                                sim_z(n,t+1)= ip
                                sim_z_preinnov(n,t)= i
                                sim_z_postinnov(n,t)= ip
                                sim_was_acquired_dummy(n,t)= 0
                                sim_rec_offer_dummy(n,t)= 1
                                sim_has_acquired_dummy(n,t)= 0
                                sim_equity_ratio(n,t)=equity_ratio(iip,jj,i,j)


                            endif
                        else ! if we receive no offer

                            sim_owner(n,t+1)=sim_owner(n,t)
                            sim_z(n,t+1)= ip
                            sim_z_preinnov(n,t)= i
                            sim_z_postinnov(n,t)= ip
                            sim_was_acquired_dummy(n,t)= 0
                            sim_rec_offer_dummy(n,t)= 0
                            sim_has_acquired_dummy(n,t)= 0
                            sim_equity_ratio(n,t)=-1d0


                        endif
                        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    else ! is target and perfect info
                        if (offermade_perfect(iip,jj,ip,j) == 1) then ! if the other firm made us an offer
                            if (gacp_eqm_perfect(iip,jj,ip,j) == 1) then  ! if we accept the offer

                                sim_owner(n,t+1)=sim_owner(n,t)+1
                                sim_z(n,t+1)=1
                                sim_z_preinnov(n,t)= i
                                sim_z_postinnov(n,t)= ip
                                sim_was_acquired_dummy(n,t)= 1
                                sim_rec_offer_dummy(n,t)= 1
                                sim_has_acquired_dummy(n,t)= 0
                                sim_equity_ratio(n,t)=equity_ratio_perfect(iip,jj,ip,j)


                            else ! if we reject the offer

                                sim_owner(n,t+1)=sim_owner(n,t)
                                sim_z(n,t+1)= ip
                                sim_z_preinnov(n,t)= i
                                sim_z_postinnov(n,t)= ip
                                sim_was_acquired_dummy(n,t)= 0
                                sim_rec_offer_dummy(n,t)= 1
                                sim_has_acquired_dummy(n,t)= 0
                                sim_equity_ratio(n,t)=equity_ratio_perfect(iip,jj,ip,j)

                            endif
                        else  !! if we receive no offer

                            sim_owner(n,t+1)=sim_owner(n,t)
                            sim_z(n,t+1)= ip
                            sim_z_preinnov(n,t)= i
                            sim_z_postinnov(n,t)= ip
                            sim_was_acquired_dummy(n,t)= 0
                            sim_rec_offer_dummy(n,t)= 0
                            sim_has_acquired_dummy(n,t)= 0
                            sim_equity_ratio(n,t)=-1d0


                        endif
                        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    endif

                endif

            endif
        enddo
    enddo

    !time_dummy(nsim1*nt1,nt1)


    deallocate(initial_z_shock)
    deallocate(psi_shock)
    deallocate(eta_shock)
    deallocate(innov_shock)
    deallocate(innov_other_shock)
    deallocate(z_other_shock)
    deallocate(tiebreaking_shock)
    deallocate(omega_shock)
    deallocate(merger_interp_shock)


    call cpu_time(finish)
    print '("Simulation Time = ",f6.3," seconds.")',finish-start





    allocate(innov(nsim1,nt1))
    allocate(innov_back_dp(nsim1,nt1))
    allocate(output_for_stata(nsim1*nt1,9))
    allocate(output_for_stata_real(nsim1*nt1,1))
    allocate(innov_stock(nsim1*nt1,1))
    allocate(innov_stock2(nsim1*nt1,1))
    allocate(innov_stock3(nsim1*nt1,1))
    allocate(innov_back_dp_vector(nsim1*nt1,2))


    mu_lowerbar = -mu_ub/2d0


    pat_dep= 16d0/17d0

    innov=0d0



    do n = 1,nsim1
        do t  = 1,nt1

            if (sim_z_postinnov(n,t)/=-1) then
                innov(n,t) = log(za(sim_z_postinnov(n,t))) - rhoa2 *  log(za(sim_z_preinnov(n,t))) - mu_lowerbar
            end if

        end do
    end do



    innov_back_dp=innov

    do n = 1,nsim1
        do t  = 2,nt1

            if (sim_owner(n,t)==sim_owner(n,t-1)  .and.sim_z_postinnov(n,t)/=-1) then
                innov_back_dp(n,t) = innov_back_dp(n,t-1)*pat_dep+innov(n,t)
            end if

        end do
    end do


    output_for_stata=0
    output_for_stata_real=0d0

    innov_back_dp_vector=0d0

    obs = 1
    do n = 1,nsim1
        do t  = 1,nt1
            output_for_stata(obs,1)=t
            output_for_stata(obs,2)=n
            output_for_stata(obs,3)=sim_owner(n,t)
            output_for_stata(obs,4)=sim_z(n,t)
            output_for_stata(obs,5)=sim_z_preinnov(n,t)
            output_for_stata(obs,6)=sim_z_postinnov(n,t)
            output_for_stata(obs,7)=sim_was_acquired_dummy(n,t)
            output_for_stata(obs,8)=sim_rec_offer_dummy(n,t)
            output_for_stata(obs,9)=sim_has_acquired_dummy(n,t)
            output_for_stata_real(obs,1)=sim_equity_ratio(n,t)
            obs = obs + 1
        end do
    end do







    if (estimate==0.and.comparastat==0) then

        OPEN(1,file='output_for_stata.txt') ; WRITE(1,"(I10)") output_for_stata


        OPEN(1,file='output_for_stata_real.txt') ; WRITE(1,"(f20.8)") output_for_stata_real
    end if




    obs=1

    !reshape to long vector and figure out the lagged stock value
    do n = 1,nsim1
        do t  = 1,nt1
            innov_back_dp_vector(obs,1)=innov_back_dp(n,t)

            if (t/=1) then
                if (sim_owner(n,t)==sim_owner(n,t-1)) then 
                    innov_back_dp_vector(obs,2)=innov_back_dp(n,t-1)
                end if
            end if
            obs = obs + 1
        end do
    end do

    if (estimate==0.and.comparastat==0) then

        OPEN(1,file='innov_back_dp_vector.txt') ; WRITE(1,"(f20.8)") innov_back_dp_vector
        OPEN(1,file='innov_back_dp.txt') ; WRITE(1,"(f20.8)") innov_back_dp
    end if
    call cpu_time(start)

    regdim=count(innov_back_dp_vector(:,1)/=0d0)
    innov_stock3=0d0
    innov_stock3(1:regdim,1)=pack(innov_back_dp_vector(:,1), innov_back_dp_vector(:,1)/= 0d0)
    mean_log_patstock=sum(innov_stock3(1:regdim,1))/regdim
    std_log_patstock=sqrt(sum((innov_stock3(1:regdim,1)-mean_log_patstock)**2d0)/regdim)


    regdim=count(innov_back_dp_vector(:,2)/=0d0)

    innov_stock=0d0
    innov_stock(1:regdim,1)=pack(innov_back_dp_vector(:,2), innov_back_dp_vector(:,2)/= 0d0) 



    allocate(norm_log_patstock(regdim,1))
    allocate(norm_log_patstock_sq(regdim,1))
    allocate(acq_dummy(regdim,1))
    acq_dummy=0d0
    acq_dummy(:,1)=real(pack(output_for_stata(:,7), innov_back_dp_vector(:,2)/= 0d0))


    norm_log_patstock(:,1)=(innov_stock(1:regdim,1)-mean_log_patstock)/std_log_patstock*2.098437+4.792221

    norm_log_patstock_sq=norm_log_patstock**2d0



    where (acq_dummy<0d0)
        acq_dummy=0d0    
    end where



    allocate(Xmatrix(regdim,3))
    allocate(Ymatrix(regdim,1))
    allocate(y_vector(regdim,1))
    allocate(XWX(3,3))
    allocate(XWXinv(3,3))
    allocate(XWY(3,1))
    allocate(bbeta_long(3,1))


    y_vector(:,1)=acq_dummy(:,1)

    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=norm_log_patstock(:,1)
    Xmatrix(:,3)=norm_log_patstock_sq(:,1)

    alpha_mul=1d0
    beta_mul=0d0

    CALL DGEMM('T','N',3,3,regdim,alpha_mul,Xmatrix,regdim,Xmatrix,regdim,beta_mul,XWX,3)


    call findinv(XWX, XWXinv, 3, errorflag)



    CALL DGEMM('T','N',3,1,regdim,alpha_mul,Xmatrix,regdim,y_vector,regdim,beta_mul,XWY,3)


    bbeta_long=matmul(XWXinv,XWY)



    !pause
    deallocate(Xmatrix)
    deallocate(Ymatrix)
    deallocate(y_vector)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(XWY)
    deallocate(norm_log_patstock)
    deallocate(norm_log_patstock_sq)
    deallocate(acq_dummy)

    !pause
    regdim=count(innov_back_dp_vector(:,1)/=0d0.and.output_for_stata(:,8)==1)
    innov_stock3=0d0
    innov_stock3(1:regdim,1)=pack(innov_back_dp_vector(:,1), innov_back_dp_vector(:,1)/= 0d0.and.output_for_stata(:,8)==1)
    mean_log_patstock=sum(innov_stock3(1:regdim,1))/regdim
    std_log_patstock=sqrt(sum((innov_stock3(1:regdim,1)-mean_log_patstock)**2d0)/regdim)

    regdim=count(innov_back_dp_vector(:,2)/=0d0.and.output_for_stata(:,8)==1)
    innov_stock2=0d0


    innov_stock2(1:regdim,1)=pack(innov_back_dp_vector(:,2), innov_back_dp_vector(:,2)/= 0d0.and.output_for_stata(:,8)==1)


    allocate(norm_log_patstock(regdim,1))
    allocate(acq_dummy(regdim,1))
    acq_dummy=0d0
    acq_dummy(:,1)=real(pack(output_for_stata(:,7), innov_back_dp_vector(:,2)/= 0d0.and.output_for_stata(:,8)==1))


    norm_log_patstock(:,1)=(innov_stock2(1:regdim,1)-mean_log_patstock)/std_log_patstock*1.901221+4.995369


    where (acq_dummy<0d0)
        acq_dummy=0d0    
    end where



    allocate(Xmatrix(regdim,2))
    allocate(Ymatrix(regdim,1))
    allocate(y_vector(regdim,1))
    allocate(XWX(2,2))
    allocate(XWXinv(2,2))
    allocate(XWY(2,1))
    allocate(bbeta_long3(2,1))


    y_vector(:,1)=acq_dummy(:,1)

    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=norm_log_patstock(:,1)

    CALL DGEMM('T','N',2,2,regdim,alpha_mul,Xmatrix,regdim,Xmatrix,regdim,beta_mul,XWX,2)


    call findinv(XWX, XWXinv, 2, errorflag)


    CALL DGEMM('T','N',2,1,regdim,alpha_mul,Xmatrix,regdim,y_vector,regdim,beta_mul,XWY,2)


    bbeta_long3=matmul(XWXinv,XWY)


    deallocate(Xmatrix)
    deallocate(Ymatrix)
    deallocate(y_vector)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(XWY)
    deallocate(norm_log_patstock)
    deallocate(acq_dummy)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Equity Ratio Simulation Regression
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    regdim=count(innov_back_dp_vector(:,1)/=0d0.and.output_for_stata(:,8)==1)
    innov_stock3=0d0
    innov_stock3(1:regdim,1)=pack(innov_back_dp_vector(:,1), innov_back_dp_vector(:,1)/= 0d0.and.output_for_stata(:,8)==1)
    mean_log_patstock=sum(innov_stock3(1:regdim,1))/regdim
    std_log_patstock=sqrt(sum((innov_stock3(1:regdim,1)-mean_log_patstock)**2d0)/regdim)

    regdim=count(innov_back_dp_vector(:,2)/=0d0.and.output_for_stata(:,8)==1)
    innov_stock2=0d0

    innov_stock2(1:regdim,1)=pack(innov_back_dp_vector(:,2), innov_back_dp_vector(:,2)/= 0d0.and.output_for_stata(:,8)==1)



    allocate(norm_log_patstock(regdim,1))
    allocate(acq_dummy(regdim,1))


    acq_dummy=0d0  !here we input equity ratio value into acq_dummy
    acq_dummy(:,1)=pack(output_for_stata_real(:,1), innov_back_dp_vector(:,2)/= 0d0.and.output_for_stata(:,8)==1)



    norm_log_patstock(:,1)=(innov_stock2(1:regdim,1)-mean_log_patstock)/std_log_patstock*1.924635 +5.028083


    allocate(Xmatrix(regdim,2))
    allocate(Ymatrix(regdim,1))
    allocate(y_vector(regdim,1))
    allocate(XWX(2,2))
    allocate(XWXinv(2,2))
    allocate(XWY(2,1))
    allocate(bbeta_long2(2,1))


    y_vector(:,1)=acq_dummy(:,1)

    Xmatrix(:,1)=1d0
    Xmatrix(:,2)=norm_log_patstock(:,1)

    CALL DGEMM('T','N',2,2,regdim,alpha_mul,Xmatrix,regdim,Xmatrix,regdim,beta_mul,XWX,2)


    call findinv(XWX, XWXinv, 2, errorflag)


    CALL DGEMM('T','N',2,1,regdim,alpha_mul,Xmatrix,regdim,y_vector,regdim,beta_mul,XWY,2)


    bbeta_long2=matmul(XWXinv,XWY)



    deallocate(Xmatrix)
    deallocate(Ymatrix)
    deallocate(y_vector)
    deallocate(XWX)
    deallocate(XWXinv)
    deallocate(XWY)




    deallocate(norm_log_patstock)
    deallocate(acq_dummy)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!




    deallocate(sim_owner)
    deallocate(sim_z)
    deallocate(sim_z_preinnov)
    deallocate(sim_z_postinnov)
    deallocate(sim_was_acquired_dummy)
    deallocate(sim_rec_offer_dummy)
    deallocate(sim_has_acquired_dummy)
    deallocate(sim_equity_ratio)







    deallocate(innov)
    deallocate(innov_back_dp)
    deallocate(output_for_stata)
    deallocate(output_for_stata_real)
    deallocate(innov_stock)
    deallocate(innov_stock2)
    deallocate(innov_stock3)
    deallocate(innov_back_dp_vector)






    modelmom2(1,1)=mean_equity_ratio
    modelmom2(2,1)=bbeta_long2(2,1)
    modelmom2(3,1)=bbeta_long(2,1)
    modelmom2(4,1)=bbeta_long(3,1)
    modelmom2(5,1)=avg_merger_prob/2d0
    modelmom2(6,1)=bbeta_long3(2,1) 
    modelmom2(7,1)=avg_realized_gain-1d0
    modelmom2(8,1)=avg_value_loss
    modelmom2(9,1) = avg_rel_size
    modelmom2(10,1) = avg_rnd_intensity
    modelmom2(11,1) = bbeta7(2,1)
    modelmom2(12,1) = std_log_mve/mean_log_mve
    modelmom2(13,1)=avg_growth_rate
    modelmom2(14,1)=psi+ avg_merger_prob/2d0



   


    if (estimate==0 .and. comparastat==0) then
        OPEN(1,file='merger_prob.txt') ; WRITE(1,"(f20.8)") merger_prob
        OPEN(1,file='conditional_accept.txt') ; WRITE(1,"(f20.8)") conditional_accept
        OPEN(1,file='equity_ratio_2.txt') ; WRITE(1,"(f20.8)") equity_ratio_2
        OPEN(1,file='equity_ratio.txt') ; WRITE(1,"(f20.8)") equity_ratio
        OPEN(1,file='equity_ratio_perfect.txt') ; WRITE(1,"(f20.8)") equity_ratio_perfect
        OPEN(1,file='corner_cash.txt') ; WRITE(1,"(f20.8)") corner_cash
        OPEN(1,file='corner_equity.txt') ; WRITE(1,"(f20.8)") corner_equity
        OPEN(1,file='interior.txt') ; WRITE(1,"(f20.8)") interior
        OPEN(1,file='new_conditional_accept.txt') ; WRITE(1,"(f20.8)") new_conditional_accept
        OPEN(1,file='new_merger_prob.txt') ; WRITE(1,"(f20.8)") new_merger_prob
        OPEN(1,file='new_merger_prob_unconditional.txt') ; WRITE(1,"(f20.8)") new_merger_prob_unconditional
        OPEN(1,file='new_offer_prob.txt') ; WRITE(1,"(f20.8)") new_offer_prob
        OPEN(1,file='new_corner_equity.txt') ; WRITE(1,"(f20.8)") new_corner_equity
        OPEN(1,file='new_corner_cash.txt') ; WRITE(1,"(f20.8)") new_corner_cash
        OPEN(1,file='new_interior.txt') ; WRITE(1,"(f20.8)") new_interior
        OPEN(1,file='new_nonoffer.txt') ; WRITE(1,"(f20.8)") new_nonoffer
        OPEN(1,file='new_equity_ratio_2.txt') ; WRITE(1,"(f20.8)") new_equity_ratio_2
        OPEN(1,file='Fa.txt') ; WRITE(1,"(f20.8)")  Fa
        OPEN(1,file='Ft.txt') ; WRITE(1,"(f20.8)")  Ft
        OPEN(1,file='new_equity_ratio_1.txt') ; WRITE(1,"(f20.8)") new_equity_ratio_1
        OPEN(1,file='realized_gain_det.txt') ; WRITE(1,"(f20.8)") realized_gain_det
        OPEN(1,file='realized_gain_det_helper.txt') ; WRITE(1,"(f20.8)") realized_gain_det_helper
        OPEN(1,file='equity_ratio_acq.txt') ; WRITE(1,"(f20.8)") equity_ratio_acq
        OPEN(1,file='normalized_zt.txt') ; WRITE(1,"(f20.8)") normalized_zt
        OPEN(1,file='normalized_z.txt') ; WRITE(1,"(f20.8)") normalized_z
        OPEN(1,file='normalized_zt2.txt') ; WRITE(1,"(f20.8)") normalized_zt2
        OPEN(1,file='v_prod.txt') ; WRITE(1,"(f20.8)") v_prod
        OPEN(1,file='fixedcost_ratio.txt') ; WRITE(1,"(f20.8)") fixedcost_ratio
        OPEN(1,file='new_realized_gain.txt') ; WRITE(1,"(f20.8)") new_realized_gain
        OPEN(1,file='realized_gain_helper.txt') ; WRITE(1,"(f20.8)") realized_gain_helper
        OPEN(1,file='new_realized_gain_dollar.txt') ; WRITE(1,"(f20.8)") new_realized_gain_dollar

        OPEN(1,file='tar_ann_ret.txt') ; WRITE(1,"(f20.8)") tar_ann_ret
        OPEN(1,file='acq_ann_ret.txt') ; WRITE(1,"(f20.8)") acq_ann_ret
        OPEN(1,file='F_acq.txt') ; WRITE(1,"(f20.8)") F_acq
        OPEN(1,file='F_tar.txt') ; WRITE(1,"(f20.8)") F_tar
        OPEN(1,file='growth_store.txt') ; WRITE(1,"(f20.8)") growth_store


        OPEN(1,file='modelmom2.txt') ;         WRITE(1,"(f20.8)") modelmom2




    end if


    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  DECOMPOSITION
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!




    if (decomposition==1) then

        do decomp_case=1,6


            decomp_diffv = 1000d0
            decomp_counter = 1

            do while (decomp_diffv>tolerv.and.decomp_counter<maxiter2)


                !$OMP PARALLEL PRIVATE(i,j,jj,gap,ttt,i_lb,i_ub,p_lb,p_ub)
                !$OMP DO collapse(3) SCHEDULE(dynamic)
                do i=1,numz
                    do j=1,numz
                        do jj=1,numx


                            gap =minval(abs(f(i,j)-za1(:,1)),1)
                            ttt =minloc(abs(f(i,j)-za1(:,1)),1)
                            if (f(i,j)>za(ttt)) then
                                i_lb=ttt
                                i_ub=min(ttt+1,numz)
                                if (i_ub>i_lb) then
                                    p_lb=1d0-gap/(za(i_ub)-za(i_lb))
                                    p_ub=1d0-p_lb
                                else
                                    p_lb=1d0
                                    p_ub=0d0
                                end if
                            else
                                i_lb=max(ttt-1,1)
                                i_ub=ttt
                                if (i_ub>i_lb) then
                                    p_lb=gap/(za(i_ub)-za(i_lb))
                                    p_ub=1d0-p_lb
                                else
                                    p_lb=0d0
                                    p_ub=1d0
                                end if
                            end if
                            V_interp(i,j,jj)=V(i_lb,jj)*p_lb+V(i_ub,jj)*p_ub


                        end do
                    end do
                end do
                !$omp end  do
                !$omp end parallel


                if (decomp_case == 2 .or. decomp_case == 5 .or. decomp_case == 6) then
                    !acquirer problem

                    
                    !$OMP PARALLEL PRIVATE(i,j,ii,jj,inn,iii,ss,cc,cc1,cc2,cc3,g1,g2,temp1,temp2,temp1b,temp2b,vv_h,vv_m,vv_l,vv,vvc,vvs,idxc,idxs,p2,x_vp)
                    !$OMP DO collapse(4) SCHEDULE(dynamic)
                    do i=1,numz  !acquirer
                        do j=1,numx !acquirer
                            do ii=1,numz !observed zt_hat
                                do jj=1,numx

                                    ss=reshape(s1,(/3,nums/),s1,(/2,1/))

                                    !print *,'ss(1,:)',ss(1,:)
                                    !  pause
                                    cc=0d0 !zeros(3,nums)

                                    cc(1,:) = c_threshold2(i,j,min(ii+1,numz),jj,:)
                                    cc(2,:) = c_threshold2(i,j,ii,jj,:)
                                    cc(3,:) = c_threshold2(i,j,max(ii-1,1),jj,:)
                                    where (cc<0d0)
                                        cc=0d0
                                    end where


                                    cc1(1,:) = cc(1,:)
                                    cc1(2,:) = cc(1,:)
                                    cc1(3,:) = cc(1,:)

                                    cc2(1,:) = cc(2,:)
                                    cc2(2,:) = cc(2,:)
                                    cc2(3,:) = cc(2,:)

                                    cc3(1,:) = cc(3,:)
                                    cc3(2,:) = cc(3,:)
                                    cc3(3,:) = cc(3,:)

                                    where (cc>=cc1)
                                        gacp_h(i,j,ii,jj,:,:)=1d0  !check this
                                    elsewhere
                                        gacp_h(i,j,ii,jj,:,:)=0d0  !check this
                                    end where

                                    where (cc>=cc2)
                                        gacp_m(i,j,ii,jj,:,:)=1d0
                                    elsewhere
                                        gacp_m(i,j,ii,jj,:,:)=0d0  !check this

                                    end where

                                    where (cc>=cc3)
                                        gacp_l(i,j,ii,jj,:,:)=1d0
                                    elsewhere
                                        gacp_l(i,j,ii,jj,:,:)=0d0  !check this

                                    end where

                                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                                    ! Imperfect Information Case
                                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                                    iii = min(ii+1,numz)
                                    x_vp = ((zt(iii)/zt(i))**phi * (zt(iii)+zt(i))**phi2)
                                    g1=LAMBDA1*( 1d0-cc/(cc+ss*(  max(0d0,pi*f(i,iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)*fv)  )     ) )**gammae2*f(i,iii) 

                                    where (g1<0d0)
                                        g1=0d0
                                    end where



                                    g2=LAMBDA1*( 1d0-cc/(cc+ss*(   max(0d0,pi*f(max(i-1,1),iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)*fv)  )       ) )**gammae2*f(max(i-1,1),iii) 
                                    where (g2<0d0)
                                        g2=0d0
                                    end where

                                    temp1=pi*f(i,iii)-x(j)*x_vp-g1
                                    temp2=pi*f(max(i-1,1),iii)-x(j)*x_vp-g2
                                    temp1b=(-cc+(1d0-ss)*(temp1+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)))
                                    temp2b=(-cc+(1d0-ss)*(temp2+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)))


                                    where (temp1b<1e-10)
                                        temp1b=1e-10
                                    end where
                                    where (temp2b<1e-10)
                                        temp2b=1e-10
                                    end where


                                    vv_h=gacp_h(i,j,ii,jj,:,:)*(kesi*temp1b(:,:)**(1d0-omega)/(1d0-omega)+(1d0-kesi)*temp2b(:,:)**(1d0-omega)/(1d0-omega))&
                                        +(1d0- gacp_h(i,j,ii,jj,:,:))*((pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))




                                    iii=ii
                                    x_vp = ((zt(iii)/zt(i))**phi * (zt(iii)+zt(i))**phi2)
                                    g1=LAMBDA1*( 1d0-cc/(cc+ss*(  max(0d0,pi*f(i,iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)*fv)  )     ) )**gammae2*f(i,iii) 
                                    where (g1<0d0)
                                        g1=0d0
                                    end where
                                    g2=LAMBDA1*( 1d0-cc/(cc+ss*(   max(0d0,pi*f(max(i-1,1),iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)*fv)  )       ) )**gammae2*f(max(i-1,1),iii) 
                                    where (g2<0d0)
                                        g2=0d0
                                    end where
                                    temp1=pi*f(i,iii)-x(j)*x_vp-g1
                                    temp2=pi*f(max(i-1,1),iii)-x(j)*x_vp-g2
                                    temp1b=(-cc+(1d0-ss)*(temp1+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)))
                                    temp2b=(-cc+(1d0-ss)*(temp2+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)))

                                    where (temp1b<1e-10)
                                        temp1b=1e-10
                                    end where
                                    where (temp2b<1e-10)
                                        temp2b=1e-10
                                    end where

                                    vv_m=gacp_m(i,j,ii,jj,:,:)*(kesi*temp1b(:,:)**(1d0-omega)/(1d0-omega)+(1d0-kesi)*temp2b(:,:)**(1d0-omega)/(1d0-omega))&
                                        +(1d0- gacp_m(i,j,ii,jj,:,:))*((pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))



                                    iii = max(ii-1,1)
                                    x_vp = ((zt(iii)/zt(i))**phi * (zt(iii)+zt(i))**phi2)
                                    g1=LAMBDA1*( 1d0-cc/(cc+ss*(  max(0d0,pi*f(i,iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)*fv)  )     ) )**gammae2*f(i,iii) 
                                    where (g1<0d0)
                                        g1=0d0
                                    end where
                                    g2=LAMBDA1*( 1d0-cc/(cc+ss*(   max(0d0,pi*f(max(i-1,1),iii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)*fv)  )       ) )**gammae2*f(max(i-1,1),iii) 
                                    where (g2<0d0)
                                        g2=0d0
                                    end where
                                    temp1=pi*f(i,iii)-x(j)*x_vp-g1
                                    temp2=pi*f(max(i-1,1),iii)-x(j)*x_vp-g2
                                    temp1b=(-cc+(1d0-ss)*(temp1+(1d0-psi)/(1d0+r)*V_interp(i,iii,j)))
                                    temp2b=(-cc+(1d0-ss)*(temp2+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),iii,j)))

                                    where (temp1b<1e-10)
                                        temp1b=1e-10
                                    end where
                                    where (temp2b<1e-10)
                                        temp2b=1e-10
                                    end where

                                    vv_l=gacp_l(i,j,ii,jj,:,:)*(kesi*temp1b(:,:)**(1d0-omega)/(1d0-omega)+(1d0-kesi)*temp2b(:,:)**(1d0-omega)/(1d0-omega))&
                                        +(1d0-gacp_l(i,j,ii,jj,:,:))*((pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega))
                                    !vv_record_big_3(i,j,ii,jj,:,:,1)=vv_h
                                    !vv_record_big_3(i,j,ii,jj,:,:,2)=vv_m
                                    !vv_record_big_3(i,j,ii,jj,:,:,3)=vv_l

                                    do inn=1,numi
                                        p2(:,:)=p3(:,:,inn)
                                        !print *, p2(1,ii),p2(2,ii),p2(3,ii)
                                        !pause
                                        vv=p2(1,ii)*vv_h+p2(2,ii)*vv_m+p2(3,ii)*vv_l



                                        vvs = vv(indexc(i,j,ii,jj),indexs(i,j,ii,jj))


                                        if (offermade(i,j,ii,jj)==1) then
                                            Vanew(i,j,ii,jj)=vvs
                                        else
                                            Vanew(i,j,ii,jj)=(pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega)
                                        endif


                                        Vanew_i(i,j,ii,jj,inn)=Vanew(i,j,ii,jj)
                                    enddo

                                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                                    ! Perfect Information Case
                                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


                                    vvs = vv_m(indexc_perfect(i,j,ii,jj),indexs_perfect(i,j,ii,jj))


                                    if (offermade_perfect(i,j,ii,jj)==1) then
                                        Vanew_perfect(i,j,ii,jj) = vvs
                                    else
                                        Vanew_perfect(i,j,ii,jj) = (pi*za(i)+(1d0-psi)/(1d0+r)*V(i,j))**(1d0-omega)/(1d0-omega)
                                    endif



                                end do
                            end do
                        end do
                    end do

                    !$omp end  do
                    !$omp end parallel

                endif

                if (decomp_case == 1 .or. decomp_case == 4 .or. decomp_case == 6) then

                    !$OMP PARALLEL PRIVATE(i,j,ii,jj,inn,iii,cc,ss,c_value,s_value,temp2s,temp2new,g1s,g2s,temp1h,temp1l,temp1s,temp1_h,temp1_m,temp1_l,x_vp)
                    !$OMP DO  SCHEDULE(dynamic)
                    do i=1,numz  !acquirer
                        do j=1,numx !acquirer
                            do ii=1,numz !target
                                do jj=1,numx

                                    ss=reshape(s1,(/3,nums/),s1,(/2,1/))

                                    !print *,'ss(1,:)',ss(1,:)
                                    !  pause
                                    cc=0d0 !zeros(3,nums)

                                    cc(1,:) = c_threshold2(i,j,min(ii+1,numz),jj,:)
                                    cc(2,:) = c_threshold2(i,j,ii,jj,:)
                                    cc(3,:) = c_threshold2(i,j,max(ii-1,1),jj,:)
                                    where (cc<0d0)
                                        cc=0d0
                                    end where

                                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                                    ! Imperfect Information Case
                                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                                    do inn=1,numi


                                        iii=min(ii+1,numz)

                                        c_value=gc_i(i,j,iii,jj,inn)
                                        s_value=gs_i(i,j,iii,jj,inn)


                                        !three cases for target


                                        temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                        temp2s=max(temp2s,1e-10)  
                                        temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                        x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                        g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                        g1s=max(g1s,0d0)
                                        g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii) 
                                        g2s=max(g2s,0d0)
                                        temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)) 
                                        temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)) 
                                        temp1h = max(temp1h,1e-10)
                                        temp1l = max(temp1l,1e-10)
                                        temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                        if (offermade_i(i,j,iii,jj,inn)==1) then

                                            temp1_h =  gacp_h(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn))*temp1s+(1d0- gacp_h(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn)))*temp2new
                                        else
                                            temp1_h=temp2new
                                        end if



                                        iii=ii
                                        c_value=gc_i(i,j,iii,jj,inn)
                                        s_value=gs_i(i,j,iii,jj,inn)



                                        temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                        temp2s=max(temp2s,1e-10)
                                        temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                        x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                        g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                        g1s=max(g1s,0d0)
                                        g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii) 
                                        g2s=max(g2s,0d0)
                                        temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)) 
                                        temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)) 
                                        temp1h = max(temp1h,1e-10)
                                        temp1l = max(temp1l,1e-10)
                                        temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                        if (offermade_i(i,j,iii,jj,inn)==1) then

                                            temp1_m =  gacp_m(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn))*temp1s+(1d0- gacp_m(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn)))*temp2new
                                        else
                                            temp1_m=temp2new
                                        end if



                                        iii=max(ii-1,1)
                                        c_value=gc_i(i,j,iii,jj,inn)
                                        s_value=gs_i(i,j,iii,jj,inn)



                                        temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                        temp2s=max(temp2s,1e-10)
                                        temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                        x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                        g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                        g1s=max(g1s,0d0)
                                        g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii) 
                                        g2s=max(g2s,0d0)
                                        temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)) 
                                        temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)) 
                                        temp1h = max(temp1h,1e-10)
                                        temp1l = max(temp1l,1e-10)
                                        temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                        if (offermade_i(i,j,iii,jj,inn)==1) then

                                            temp1_l =  gacp_l(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn))*temp1s+(1d0- gacp_l(i,j,iii,jj,indexc_i(i,j,iii,jj,inn),indexs_i(i,j,iii,jj,inn)))*temp2new
                                        else
                                            temp1_l=temp2new
                                        end if



                                        Vtnew_h(i,j,ii,jj)=temp1_h
                                        Vtnew_m(i,j,ii,jj)=temp1_m
                                        Vtnew_l(i,j,ii,jj)=temp1_l


                                        Vtnew_h_i(i,j,ii,jj,inn)=temp1_h
                                        Vtnew_m_i(i,j,ii,jj,inn)=temp1_m
                                        Vtnew_l_i(i,j,ii,jj,inn)=temp1_l

                                    enddo


                                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                                    ! Perfect Information Case
                                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                                    iii=ii
                                    c_value=cc(indexc_perfect(i,j,iii,jj),indexs_perfect(i,j,iii,jj))
                                    s_value=ss(indexc_perfect(i,j,iii,jj),indexs_perfect(i,j,iii,jj))


                                    temp2s=outside_option*(pi*zt(ii)+(1d0-psi)/(1d0+r)*V(ii,jj))
                                    temp2s=max(temp2s,1e-10)
                                    temp2new=temp2s**(1d0-omega2)/(1d0-omega2)

                                    x_vp = ((zt(ii)/zt(i))**phi * (zt(ii)+zt(i))**phi2)
                                    g1s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(   max(0d0,pi*f(i,ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)*fv)      )) )**gammae2*f(i,ii) 
                                    g1s=max(g1s,0d0)
                                    g2s=LAMBDA1*( 1d0-c_value/(c_value+s_value*(     max(0d0,pi*f(max(i-1,1),ii)-x(j)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)*fv)       )) )**gammae2*f(max(i-1,1),ii) 
                                    g2s=max(g2s,0d0)
                                    temp1h = c_value+s_value*(pi*f(i,ii)-x(j)*x_vp-g1s+(1d0-psi)/(1d0+r)*V_interp(i,ii,j)) 
                                    temp1l = c_value+s_value*(pi*f(max(i-1,1),ii)-x(j)*x_vp-g2s+(1d0-psi)/(1d0+r)*V_interp(max(i-1,1),ii,j)) 
                                    temp1h = max(temp1h,1e-10)
                                    temp1l = max(temp1l,1e-10)
                                    temp1s = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)

                                    if (offermade_perfect(i,j,iii,jj)==1) then

                                        temp1_m =  gacp_m(i,j,iii,jj,indexc_perfect(i,j,iii,jj),indexs_perfect(i,j,iii,jj))*temp1s+(1d0- gacp_m(i,j,iii,jj,indexc_perfect(i,j,iii,jj),indexs_perfect(i,j,iii,jj)))*temp2new
                                    else
                                        temp1_m=temp2new
                                    end if

                                    Vtnew_perfect(i,j,ii,jj) = temp1_m

                                end do
                            end do
                        end do
                    end do
                    !
                    !            !
                    !$omp end  do
                    !$omp end parallel

                endif


                do i=1,numz  ! my pre-innovation productivity
                    do j=1,numx ! my merger cost

                        do inn=1,numi ! my innovation choice

                            do ii=1,numz ! other pre-innovation productivity
                                do jj=1,numx ! other merger cost

                                    if (decomp_case == 1 .or. decomp_case == 4 .or. decomp_case == 6) then
                                        tempt = 0d0
                                        do k = -1,1 ! my innovation realization
                                            do m = -1,1 ! other innovation realization
                                                ip = max(min(i+k,numz),1)
                                                iip = max(min(ii+m,numz),1)

                                                if (decomp_case == 1) then
                                                    if (k == -1) then
                                                        tempt = tempt + (info_asym*Vtnew_h_i(iip,jj,ip,j,inn) + (1d0-info_asym)*Vtnew_perfect(iip,jj,ip,j))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                                    elseif (k == 0) then
                                                        tempt = tempt + (info_asym*Vtnew_m_i(iip,jj,ip,j,inn) + (1d0-info_asym)*Vtnew_perfect(iip,jj,ip,j))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                                    else
                                                        tempt = tempt + (info_asym*Vtnew_l_i(iip,jj,ip,j,inn) + (1d0-info_asym)*Vtnew_perfect(iip,jj,ip,j))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                                    endif
                                                else
                                                    if (k == -1) then
                                                        tempt = tempt + (Vtnew_perfect(iip,jj,ip,j))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                                    elseif (k == 0) then
                                                        tempt = tempt + (Vtnew_perfect(iip,jj,ip,j))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                                    else
                                                        tempt = tempt + (Vtnew_perfect(iip,jj,ip,j))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                                    endif
                                                endif
                                            enddo
                                        enddo
                                    endif

                                    if (decomp_case == 2 .or. decomp_case == 5 .or. decomp_case == 6) then
                                        tempa = 0d0
                                        do k = -1,1 ! my innovation realization
                                            do m = -1,1 ! other innovation realization
                                                ip = max(min(i+k,numz),1)
                                                iip = max(min(ii+m,numz),1)

                                                if (decomp_case == 2) then
                                                    tempa = tempa + (info_asym*Vanew_i(ip,j,ii,jj,i_policy(ii,jj)) + (1d0-info_asym)*Vanew_perfect(ip,j,iip,jj))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                                else
                                                    tempa = tempa + (Vanew_perfect(ip,j,iip,jj))*p3(2-k,i,inn)*p3(2-m,ii,i_policy(ii,jj))
                                                endif
                                            enddo
                                        enddo
                                    endif

                                    temp_alone = 0d0
                                    do k = -1,1 ! my innovation realization
                                        ip = max(min(i+k,numz),1)
                                        temp_alone = temp_alone + (pi*za(ip)+(1d0-psi)/(1d0+r)*V(ip,j))*p3(2-k,i,inn)
                                    enddo

                                    if (decomp_case == 1 .or. decomp_case == 4) then ! target only case (unable to be an acquirer)
                                        if (i < ii) then
                                            V_meet_i(ii,jj,inn) = tempt
                                        elseif (i == ii) then
                                            V_meet_i(ii,jj,inn) = (tempt + temp_alone)/2d0
                                        else
                                            V_meet_i(ii,jj,inn) = temp_alone
                                        endif
                                    elseif (decomp_case == 2 .or. decomp_case == 5) then ! acquirer only case (unable to be a target)
                                        if (i < ii) then
                                            V_meet_i(ii,jj,inn) = temp_alone
                                        elseif (i == ii) then
                                            V_meet_i(ii,jj,inn) = (temp_alone + tempa)/2d0
                                        else
                                            V_meet_i(ii,jj,inn) = tempa
                                        endif
                                    elseif (decomp_case == 6) then
                                        if (i < ii) then
                                            V_meet_i(ii,jj,inn) = tempt
                                        elseif (i == ii) then
                                            V_meet_i(ii,jj,inn) = (tempt + tempa)/2d0
                                        else
                                            V_meet_i(ii,jj,inn) = tempa
                                        endif
                                    else ! always alone case (unable to be an acquirer or a target)
                                        V_meet_i(ii,jj,inn) = temp_alone
                                    endif

                                enddo
                            enddo



                            Vnew_i(i,j,inn) = eta*(sum(V_meet_i(:,:,inn)*Fs(:,:)/sum(Fs))) + (1d0-eta)*temp_alone - (i_cost(inn)*za(i))

                        enddo
                        Vnew(i,j) = Vnew_i(i,j,i_policy_new(i,j))
                    enddo
                enddo

                decomp_diffv=maxval(abs(Vnew-V)) 
                decomp_counter=decomp_counter+1

                Va=Vanew 
                Vt=Vtnew 
                V=Vnew*vfi_update+ V*(1d0-vfi_update)
            enddo

            print *, 'decomp_case decomp_counter decomp_diffv', decomp_case, decomp_counter, decomp_diffv

            if (decomp_case == 1) then
                OPEN(1,file='V_decomp_tar_only.txt') ; WRITE(1,"(f20.8)") V
            elseif  (decomp_case == 2) then
                OPEN(1,file='V_decomp_acq_only.txt') ; WRITE(1,"(f20.8)") V
            elseif  (decomp_case == 3) then
                OPEN(1,file='V_decomp_alone.txt') ; WRITE(1,"(f20.8)") V
            elseif  (decomp_case == 4) then
                OPEN(1,file='V_decomp_tar_no_asym.txt') ; WRITE(1,"(f20.8)") V
            elseif (decomp_case == 5) then
                OPEN(1,file='V_decomp_acq_no_asym.txt') ; WRITE(1,"(f20.8)") V
            else
                OPEN(1,file='V_decomp_all_no_asym.txt') ; WRITE(1,"(f20.8)") V
            endif

        enddo

    end if













    end subroutine makemoments



    subroutine equation2(cv,outside_value,s_idx,i_idx,ii_idx,j_idx)

    use parameters
    use globals
    use tools
    implicit none

    real*8,intent(in):: cv,outside_value
    integer::s_idx,i_idx,ii_idx,j_idx
    real*8::g1s0,g2s0


    x_vp = ((zt(ii_idx)/zt(i_idx))**phi * (zt(ii_idx)+zt(i_idx))**phi2)
    g1s0=LAMBDA1*( 1d0-cv/(cv+s(s_idx)*(  max(0d0,pi*f(i_idx,ii_idx)-x(j_idx)*x_vp+(1d0-psi)/(1d0+r)*V_interp(i_idx,ii_idx,j_idx)*fv)      )) )**gammae2*f(i_idx,ii_idx)
    g1s0=max(g1s0,0d0)
    g2s0=LAMBDA1*( 1d0-cv/(cv+s(s_idx)*(  max(0d0,pi*f(max(i_idx-1,1),ii_idx)-x(j_idx)*x_vp+(1d0-psi)/(1d0+r)*V_interp(max(i_idx-1,1),ii_idx,j_idx)*fv)     )) )**gammae2*f(max(i_idx-1,1),ii_idx)
    g2s0=max(g2s0,0d0)
    temp1h = cv+s(s_idx)*(pi*f(i_idx,ii_idx)-x(j_idx)*x_vp-g1s0+(1d0-psi)/(1d0+r)*V_interp(i_idx,ii_idx,j_idx))
    temp1l = cv+s(s_idx)*(pi*f(max(i_idx-1,1),ii_idx)-x(j_idx)*x_vp-g2s0+(1d0-psi)/(1d0+r)*V_interp(max(i_idx-1,1),ii_idx,j_idx))
    temp1h = max(temp1h,1e-10)
    temp1l = max(temp1l,1e-10)
    temp1_record = kesi*temp1h**(1d0-omega2)/(1d0-omega2) + (1d0-kesi)*temp1l**(1d0-omega2)/(1d0-omega2)


    value_diff = temp1_record - outside_value




    end subroutine equation2




    subroutine giveindex(z_idx,x_idx)

    use parameters
    use globals
    use tools
    implicit none

    integer,intent(in):: z_idx, x_idx


    zx_idx=(x_idx-1)*numz + z_idx

    end subroutine giveindex




    subroutine zinterp(za_idx,zt_idx,p_lb0, p_ub0,i_lb0,i_ub0)

    use parameters
    use globals
    use tools
    implicit none

    integer,intent(in):: za_idx, zt_idx
    integer, intent(out)::i_lb0,i_ub0
    real*8, intent(out)::p_lb0, p_ub0
    real*8::gap0
    integer::tt
    gap0=minval(abs(f(za_idx,zt_idx)-za1(:,1)),1)
    tt =minloc(abs(f(za_idx,zt_idx)-za1(:,1)),1)
    if (f(za_idx,zt_idx)>za(tt)) then
        i_lb0=tt
        i_ub0=min(tt+1,numz)
        if (i_ub0>i_lb0) then

            p_lb0=1d0-gap0/(za(i_ub0)-za(i_lb0))
            p_ub0=1d0-p_lb0
        else
            p_lb0=1d0
            p_ub0=0d0
        end if
    else
        i_lb0=max(tt-1,1)
        i_ub0=tt
        if (i_ub0>i_lb0) then
            p_lb0=gap0/(za(i_ub0)-za(i_lb0))
            p_ub0=1d0-p_lb0
        else
            p_lb0=0d0
            p_ub0=1d0
        end if
    end if




    end subroutine zinterp




























    subroutine readweightmatrix(cluster, weight)

    use parameters
    use globals
    use tools
    implicit none

    integer, intent(in) :: cluster
    real(8), intent(out) :: weight(nmom, nmom)

    integer ::  errflag
    real :: temp

    real(8) :: iweightread(nmom, nmom) !iweightread(allnmom, allnmom)

    character (len=50) :: covfil

    weight=0d0

    do i=1,nmom
        weight(i,i) = 1d0
    end do


    !if (setting==1) then
    !
    !    if (cluster == 1) then
    !        covfil = "Datfils/omega.txt"   
    !        open (unit=62,file=trim(covfil),status="old")
    !    else
    !        ! covfil = "Datfils/covmx.txt"
    !        covfil = "Datfils/Wcovmx.txt"
    !        open (unit=62,file=trim(covfil),status="old")
    !    endif
    !
    !    do ii = 1, nmom
    !        read(62,*) iweightread(:, ii)
    !    enddo
    !    close (unit=62)
    !elseif (setting==2) then
    !    if (cluster == 1) then
    !        covfil = "Datfils/omega_early.txt"  
    !        open (unit=62,file=trim(covfil),status="old")
    !    else
    !        ! covfil = "Datfils/covmx.txt"
    !        covfil = "Datfils/Wcovmx_early.txt"
    !        open (unit=62,file=trim(covfil),status="old")
    !    endif
    !
    !    do ii = 1, nmom
    !        read(62,*) iweightread(:, ii)
    !    enddo
    !    close (unit=62)
    !elseif (setting==3) then
    !    if (cluster == 1) then
    !        covfil = "Datfils/omega_late.txt"   
    !        open (unit=62,file=trim(covfil),status="old")
    !    else
    !        ! covfil = "Datfils/covmx.txt"
    !        covfil = "Datfils/Wcovmx_late.txt"
    !        open (unit=62,file=trim(covfil),status="old")
    !    endif
    !
    !    do ii = 1, nmom
    !        read(62,*) iweightread(:, ii)
    !    enddo
    !    close (unit=62)
    !end if
    !
    !!      do ii = 1, nmom
    !!         read(62,*) weight(:, ii) !"(<numinfl>F26.16)"
    !!      enddo
    !!         close (unit=62)
    !
    !
    !do ii = 1, nmom
    !    do jj = 1, nmom
    !        weight(ii, jj) = iweightread(ii, jj) !iweightread(pickmom(ii), pickmom(jj))
    !    enddo
    !enddo

    end subroutine readweightmatrix






